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PART 1: INTRODUCTION-—The role 
analog computers can play in engineer- 
ing is discussed for the practicing 
engineer. Computer components 

are identified 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


Caxcutations, previously done by hand, can now be 
completed using computers with a great reduction in cal- 
culation time and a marked increase in computation pre- 
cision, While such incentives are frequently sufficient to 
justify the widespread interest in computers, the applica- 
tion of computers to the solution of engineering problems 
has had a much more fundamental effect on chemical 
and petroleum engineering. When the speed and precision 
of modern computers are coupled with their variety, ver- 
satility, and large capacity, the reasons soon become 
apparent. 


The practicing engineer now has at his disposal com- 
putational devices which enables him to handle economi- 
cally problems of a complexity which he could not have 
handled 10 years ago. As a result, the class of problems 
which the engineer can be expected to handle has been 
greatly extended. This extension has resulted in a re- 
evaluation of what can be classed engineering calculations. 
The re-evaluation is still in progress, but it appears clear 
that the engineer of today, and certainly of tomorrow, 
will be constantly associated with problems once consid- 
ered too complex, time consuming, or precision oriented 
to be amicable of solution. Consequently, computers have 
not only proven to be useful computational tools for the 
engineer, they have also provided a means for extending 
the scope of engineering calculations and, in so doing, 
have extended the scope of engineering itself. 


The increasing complexity of engineering problems is 
recognized by universities so that engineers today are 
being routinely trained in the use of computers, as they 
in the past were trained in the use of slide rules and 
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hand calculating machines. While the degree of computer 
competence which is required of a practicing engineer 
varies with his position, some knowledge of computers 
must be considered an integral part of his already varied 
skills. 

For the practicing engineer who is directly associated 
with the programing and operation of computers, a high 
degree of competence on the computer is required. A 
common situation found in many companies is that in 
which a group of experts, highly trained in the use of 
computers, is available for consultation with the practic- 
ing engineer. In this case, the practicing engineer can 
rely on these experts for assistance in programing and 
operation, In many cases, these experts will handle all of 
the details of obtaining a solution once the engineer has 
properly defined his problem. 


If experts are available, the degree of competence 
required by an engineer to solve a specific problem on a 
specific computer is inversely related to the degree of help 
he can obtain from the experts. However, frequently the 
engineer's problems are vaguely defined. Assumptions of 
varying degree must frequently be made, the reliability 
of the data examined, and, if several computers are avail- 
able, the choice of which to use must be made. Conse- 
quently, the engineer is seldom concerned only with the 
solution of a specific problem on a specific computer, and, 
even when computer experts are available, knowledge of 
the use of computers can significantly aid the practicing 
engineer. 


A background of computer knowledge: 


© Enables the engineer to communicate more effec- 
tively with the expert. It permits him to state his problem 
more clearly to the expert and to anticipate the potential 
pitfalls which the expert may encounter during the de- 
tails of solution. 


© Makes the engineer aware of the limitations of the 
available computers. This aids the engineer in making the 
choice of computer to use to solve his current problem 
and enables him to gear his method of solution so as to 
minimize the importance of the limitations. The engineer 
will also be aware of the versatility of the chosen com- 
puter and be able to decide*the feasibility of relaxing 
some of his assumptions on an individual program. 

© Perhaps most important, stimulates the engineer to 
attack problems which he could not previously consider. 
By being aware of the speed, precision, variety, versatility, 
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and capacity of computers as well as of the techniques of 
application, the engineer will not lightly brush off a 
unique problem because of its complexity. Indeed, he will 
be stimulated to examine new approaches to problems 
and new solution techniques as they may arise in his 
day-to-day activities. 





Comparison of Computers. Modern computers may be 
classified as either digital, analog, or hybrid computer: 


The digital computer has the characteristic that the 
arithmetic operations of addition, subtraction, multiplica- 
tion and division of only discrete values can be accom- 
plished very rapidly. Digital computers can be equipped 
with a large memory capacity which permits the storage 
of and the operation of an almost limitless number of 
discrete values. Moreover, the arithmetic operations can be 
performed to a high degree of precision which is generally 
fixed by the manufacturer on any particular machine. 
Thus the digital computer is in a sense the equivalent of 
a very reliable and extremely rapid desk calculator. 

In addition to arithmetic operations the digital com- 
puter can be used to make logical decisions by comparing 
two discrete values and deciding whether one is equal to, 
less than or greater than the other. The result of the 
digital’s facilities is a most versatile computational tool 
which when coupled with the numerical techniques for 
trial-and-error, integration, etc. permits the handling of a 
very wide class of problems. However, since only discrete 
values are used and calculations are done sequentially, 
considerable time may be required for problems in which 
a large number of discrete values must be used. Program- 
ing a digital computer is a very precise operation. Detailed 
instructions, generally in a very specialized language, must 
be given for every calculation. Digital computers find 
widespread application in data reduction, statistics, solu- 
tion of algebraic equations, steady-state calculations, etc., 
where discrete values occur. Numerical technique and the 
accompanying versatility of the digital computer further 
extend its range of application. 


The analog computer, in sharp contrast to the digital 
performs computations in parallel and continuously. The 
computations which can be performed include the arith- 
metic operations as well as continuous arbitrary function 
operations, such as squaring, taking the logarithm, ex- 
ponentiation, etc. The most significant feature of the 
analog, however, is the capability of performing integra- 
tion on a continuous basis. The analog computer excels 
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TABLE 1—Anclog Computer Components 
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in the solution of differential equations regardless of 
whether they are linear or nonlinear. The solution from 
an analog computer is presented in continuous form 
which is analogous to the exact solution rather than in 
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terms of discrete values which can only approximate the 
exact solution. Furthermore, the computations can be car- 
ried out at practically any speed desired. This permits 
the rapid examination of a wide range of parameters in 
a space of time that can be significantly shorter than the 
fixed time of calculation on the digital. Analogs do not, 
however, have the capacity for memory, the capacity for 
a large number of operations, nor the degree of precision 
of the digital. 

Programing an analog does not require knowledge 
of a highly specialized language as does the digital. The 
details of programing are generally quite similar to classi- 
cal methods of solution with which the engineer is famil- 
iar. Analog computers find their greatest application in 
the solution of differential equations, such as occur in the 
dynamic analysis of processing systems, in which the 
amount of algebraic and logical operations required is 
rather limited. 


A hybrid computer is the combination of a digital and 
an analog computer in which the solution of a problem 
is shared between the two. By containing both an analog 
and a digital, the solution of a problem on the hybrid 
computer can incorporate the advantages of both. The 
digital computer portion provides a large memory capac- 
ity and permits rapid logical and algebraic operations 
while the analog permits continuous integration. The 
transfer of information from one portion to the other and 
back again extends the versatility of the hybrid beyond 
that of the analog or the digital alone. Hybrid computers 
are finding application in the solution of complex prob- 
lems where the time required for a digital solution is ex- 
cessive and where analogs do not possess the capacity or 
the precision required for a solution. 


Analog Computer Components, The term analog in 
analog computer is really a misnomer. That is, the electri- 
cal circuitry used to solve a particular problem is not an 
electrical analog in the conventional sense that electrical 
current is analogous to fluid flow, voltage is analogous to 
pressure, electrical capacitance is analogous to mass, etc. 
The electrical analog computer is actually a device which 
has been designed to perform certain mathematical opera- 
tions, such as addition, multiplication, integration, etc. on 
specified voltages. The use of an analog computer conse- 
quently involves the specification and completion of the 
mathematical operations required to solve a particular 
set of equations and is entirely independent of the fact 
that the equations may describe a flow system, a heat 
exchange problem, a mass transfer problem, or a reactor 
system. 

Since most of the analog computers being used today 
are transistorized, this discussion is restricted to this type 
of analog computer. This is not a large restriction because 
most of the older vacuum tube computers operate in the 
same fashion with the exception of a few mechanical com- 
ponents. 


A list of the components found in a modern electronic 
analog computer is given in Table 1. The first column 
in Table 1 contains a symbolic representation of each 
component, the second column contains an electrical rep- 
resentation of each component and the third column shows 
the mathematical operation that each component performs. 
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Although it is not shown in Table 1, one of the most 
important components of the computer is the constant 
voltage DC power supply. Depending upon the design of 
the computer, the power supply is usually either 100 or 
10 volts. This DC voltage is called the reference or ma- 
chine voltage and its stability will determine to a large 
degree the accuracy of the solution that one obtains. 

The first element in Table 1 is an attenuator com- 
monly called a pot—an abbreviation for potentiometer. 
Two types of pots are commonly found in an analog 
computer. One is referred to as a grounded pot and the 
other an ungrounded pot. The grounded pot enables the 
operator to obtain an output voltage which is some frac- 
tion between 0 and 1 of the input voltage. The un- 
grounded pot performs the function shown in Table 1 
and is generally used in the construction of special func- 
tion circuits, Attenuators are usually adjusted manually 
to set the value of a constant K. This value is determined 
by measuring the output voltage from an attenuator with 
the voltage measuring device built into the computer. 

The next component listed in Table 1 is the high gain 
amplifier. The function of the high gain amplifier is to 
multiply an input voltage by a large constant, generally 
of the order of 10%. The high gain amplifier is seldom 
used by itself, but usually forms an integral part of other 
common analog components, The high gain amplifier is 
generally designed so that it provides very stable operation 
if the magnitude of the output voltage does not exceed 
the computer reference voltage. If the output voltage of 
the high gain amplifier exceeds the machine voltage the 
operation of the computer will be non-linear. For this 
reason, the output of a high gain amplifier during opera- 
tion should not exceed that of the computer reference 
voltage. Most computers have a built in alarm system 
which warns the operator when the output of an ampli- 
fier exceeds machine voltage. This restriction on the out- 
put voltage holds for every component in which a high 
gain amplifier appears. 

The next component of interest in Table 1 is the sum- 
ming circuit. As its name implies, the summing circuit 
adds two or more voltages and gives the resulting sum. 
The summing circuit consists of a high gain amplifier 
with a feedback resistor and three or more input resistors. 
The input resistors and the feedback resistor, in a sum- 
ming circuit, are usually fixed not variable resistors. The 
input voltage to a summing circuit is multiplied by the 
negative ratio of the resistance of the feedback resistor 
to the resistance of the input resistor. For most computers 
the value of this ratio is usually 0.1, 1 or 10, Since a sum- 
mer contains a high gain amplifier, the output from a 
summing circuit must not exceed the reference voltage. 

One of the most important components of an analog 
computer is the integrator. The circuit for an integrator 
is essentially a high gain amplifier equipped with a capaci 
tor in the feedback line and two or more input resistors, 
The output from an integrator is the negative sum of 
the initial condition and the integral with respect to com- 
puter operating time of fixed ratios of the input voltages 
as indicated in Table 1. Since the output of the integrator 
is time varying, there is a possibility that this value may 
be greater than machine voltage at sometime during the 
problem. Therefore, special precautions must be taken tc 
insure that this does not happen. These precautions wil 
be discussed later in the series. 

The components discussed to this point are common) 











called linear components and the principle of superposi- 
tion holds. The next group of components to be discussed 
are commonly called non-linear components. 


Nonlinear mathematical operations, including multi- 
plication, division, exponentiation, etc., are achieved 
through the use of specially designed components. Gen- 
erally a particular component can be used to obtain both 
the designed operation and its inverse depending on the 
particular patching used. Thus, what is termed a multi- 
plier can be used to achieve either multiplication or divi- 
sion. The details of patching and the electronic circuitry 
varies from one type of analog to another. The usually 
acceptable programing symbols are in Table 1. 

A discussion of the internal construction of the non- 
linear components will be covered in a later part. The 
purpose here is to present some of the more common 
nonlinear components and to examine their limitations 
and the precautions which must be exercised in their use. 

When it is desired to generate a function such as 
y = x*, where x is the input and y is the output, a non- 
linear component is used to achieve this operation. Tran- 
sistorized nonlinear components can be used to produce 
an output which is a series of straight-line approximations 
as indicated by the dotted lines in Fig. 1. 

Generally, for economic reasons, sufficient line segments 
are not used to closely approximate the desired function 
over the entire input range. The best response is usually 
obtained when the input is near its maximum magnitude 
(that is, the machine reference voltage) . 


A high gain amplifier is an integral part of nonlinear 
components. Because of this, the output from a non- 
linear component can not exceed the reference voltage. 
Numerical factors are included in nonlinear components 
so that when the input is near its maximum magnitude 
the output is also. 

Because certain nonlinear operations such as squaring 
may be more easily obtained than other operations, 
mathematical relationships such as Equation (1) may 
form an integral part of a nonlinear component. Perhaps 
the most common of these is the relationship used in the 
quarter-square multiplier 


yee 





(a) 


which is used to obtain a product because squaring is 
fairly easily achieved. 

Due to the internal construction of nonlinear com- 
ponents, they will draw a varying amount of current de- 
pending on the value of the input. It is, therefore, very 
important that the output from a pot never be used as 
an input to a nonlinear device since a varying current 
produces a varying voltage distribution from the pot. In- 
puts to nonlinear components should always be outputs 
from high gain amplifiers contained in other analog com- 
ponents. 

The size of an analog computer is usually gaged by 
the number of high gain amplifiers that it contains. Most 
commercial computers normally run from ten to several 
hundred high gain amplifiers, most of these high gain 
amplifiers are normally summers. A somewhat smaller 
number are integrators, which may also be used as sum- 
mers. The remainder of the high gain amplifiers are 
usually integral parts of specialized nonlinear components. 
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Fig. 1—The nonlinear relation shown by the solid line can be 
approximated by a series of linear relations. 
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Fig. 2—An analog computer circuit is developed to compute 
the height of liquid in this tank. 
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Fig. 3—Steps are combined to give the finished circuit. 


The complexity of the problem which may be solved is 
very closely related to the number of high gain amplifiers 
available on the computer. A problem is usually solved 
on an analog computer by connecting the various com- 
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Fig. 4—Numerical values are added to the circuit. 


Fig. 5—The analog circuit 
is simplified. 





ponents on a patch panel with external leads. This process 
is called patching. 

As an example of how an analog computer circuit 
may be set up we will consider the case of the change 
in height of liquid in a tank having a cross-section area A. 

A mass balance on this system yield 








Accumulation = Input — Output (2) 
A (dh/dt) = 4, — 4 (3) 


We will assume that q; is constant and that the rate of 
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flow out of the tank is linearly related to the height of 
fluid in the tank through the flow line resistance R». 


a2 = W/R, (4) 








Then our relation becomes upon substitution and re- 
arrangement 


dh/dt = (q,/A) — (h/AR,) (5) 


The analog diagram may be constructed from the dif- 
ferential equation by applying the following steps: 

Step 1. Write the differential equation with the highest 
ordered derivative by itself on the left side of the equation. 

Step 2. Begin a diagram by showing the highest ordered 
derivative at the output of a summing amplifier. 

Step 3. Perform as many integrations as are necessary to 
form all of the variables on the right side of the equation. 

Step 4. Use the generated variables to form the original 
highest ordered derivative. 

Step 5. Insert any forcing functions. 

The sequential construction of the program is shown 
in Fig. 3. 

The pot settings can be determined when the values 
of parameters in the problem are known. Assume that 





Rynit 





10 sq. ft. 
0 cu, ft./min 
min./sq. ft. 


Pot 3 = Aynitiar/10 = 0.4 





The final diagram with pot settings is given in Fig. 4 
Although the analog diagram appears satisfactory, it may 
indeed be unsatisfactory for one or both of two reasons. 

© The output from a high gain amplifier may exceed 
the reference voltage and hence yield a faulty solution. 

The output from a high gain amplifier may be so 
small that inherent limitations of the analog results in 
a faulty solution or difficulty is encountered in monitoring 
the value. 

From a knowledge of the nature of the numerical solu- 
tion, the output, h, from amplifier 2 has an initial value 
of 4 v., a final value of 10 v., and asymptotically ap- 
proaches the final value in an exponential fashion. Hence 
for a 100 v. or a 10 v. machine, amplifier 2 will not be 
overloaded, Similarly amplifier 1 with the output dh/dt 
will not be overloaded because it has an initial value of 0.6 
and decreases montonically to 0, However, for a 100 v. 
machine, the values of h and dh/dt are quite small and 
for a 10 v. machine dh/dt is quite small. Since the cir- 
cuit uses no nonlinear components, the solution will prob- 
ably be satisfactory although monito1 ing may be difficult. 

It should be noted that, if the value of the highest de- 
rivative does not need to be available for monitoring, 
step 2 above may be eliminated. In this case, the diagram 
shown in Fig. 5 is obtained for the chosen example 
with a net reduction in the number of analog components 
required for a solution. 

In the next article of this series, techniques will be pre- 
sented which permit an analog solution to be obtained 
when the aforementioned reasons prevent a direct solu- 
tion. The process by which satisfactory computer variables 
are chosen so that the outputs from amplifiers are main 
tained near, but do not exceed the reference voltage, is 
termed amplitude scaling. 
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PART 2: AMPLITUDE SCALING—The 
outputs from high gain amplifiers 
in an analog computer should be 
maintained near the machine’s 
reference voltage, but not exceed it 


Theodore W. Cadman and Theodore G. Smitl 
University of Maryland, College Park, Md. 


For SATISFACTORY OPERATION, the output from a high 
gain amplifier contained in any analog component must 
not exceed the machine reference voltage. On the other 
hand, the outputs should be maintained near their maxi- 
mum permissible values for the most satisfactory opera- 
tion of most nonlinear components and for ease and pre- 
cision of monitoring. Thus a method is needed to scale 
an analog problem so that the outputs from high gain 





amplifiers are maintained near, but do not exceed, the 
reference voltage. This process is termed amplitude scal- 
ing. 


Amplitude Scaling. This is accomplished by expressing 
equations in a normalized form, wherein the normalized 
variables are to be the outputs from high gain amplifiers 
with magnitudes which approach, but do not exceed, the 
reference voltage. The analog solution is consequently the 
solution of the normalized form. 

The process of amplitude scaling is perhaps most easily 
achieved if the reference voltage is considered to be + 1 
unit, even though it may be + 10 volts or + 100 volts. 
Scaling accomplished using this generalization is termed 
unity scaling and has several advantages over scaling in 
which the value of the reference voltage is implicitly used. 

First, unity scaling permits an analog computer pro- 
gram to be obtained which can be used on any analog 
computer regardless of the value of the reference voltage. 
Second, the normalized variables which are outputs from 


Fig. 6—Unity scaling is applied to an engineering problem before it is solved on the analog computer. 
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TABLE 2—Scaling Ampli 





yr Outputs and Potentiometer Settings 











Output Estimated 
Component or Maximum Sees 
Value a 
Amp 1 h,* h,/h,* 
Amp 2 h,* —h,/h,* 
Amp 3 h* hi /hy* 
Amp 4 h* —h,/h,* 
Pot 1 = h,(0)/h,* 
Pot 2 = 9,/A,hy* 
Pot 3 = 1/A,R, 
Pot 4 = 
Pot 5 at 
Pot 6 = 
Pot 7 en =— 
2 AR, 
1 
Pot 8 ee ew 
a aR, 








high gain amplifiers are simply the value of the original 
variable divided by the maximum magnitude of the origi- 
nal variable. Third, and perhaps most significant, the 
process of scaling non-linear components is greatly sim- 
plified. 

For example, if x and y are the inputs to a multiplier, 
the output expressed in volts for a 10-volt machine is 
xy/10. For a 100-volt machine, the output is xy/100 ex- 
pressed in volts. But in terms of unity scaling it is xy 
regardless of the machine reference voltage. The simplifi- 
cation encountered for other nonlinear components is 
similar. 


Unity scaling can be achieved by mathematical ma- 
nipulation of the original equations before any program- 
ing is attempted. However, scaling is frequently more 
easily achieved and more clearly visualized if it is done 
after some preliminary programing. In particular, the 
following seven steps have been found to be a most con- 
venient method of achieving unity scaling. 


Step 1. Write down the differential and/or algebraic 
equations whic! lem_to ed. 

Step 2. Sketch a preliminary analog diay eglect- 
ing the prebien of alm OF scaling, which connects all of the re- 
quired analog components together. Clearly mark the 


output from each amplifier and the value of each pot as 
it appears in the problem equations. 


Step 3. Make a_list of the variables that are outputs. 
ese are the variables which must be 


from amplifiers. 
anata Sead. Also list the values of the pots. These 


may be changed as a result of the amplitude scaling 
process. 


Step 4. Estimate the maximum magnitude of each of. 
the variables which are outputs from amplifiers, 

Step 5. Form the normalized or computer variables for 
the outputs from the amplifiers by dividing the problem 
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variable by its estimated maximum magnitude. Note these 
on the list made in Step 3. 


Step 6. Insert additional pots in the amplifier input 
lines and/or modify the values of the pots which were 
needed to complete Step 2 so as to compensate for the 
change in scale from one amplifier to the next in series. 
Again modify the list completed in Step 3. 


Step 7. Sketch the completed analog diagram using the 
normalized variables and the new pot settings. Note that 
this step is easily accomplished by referring to the pre- 
liminary sketch and to the list which has been completed. 


If the programer closely follows these rules and system- 
atically lists each variable and the pot settings during the 
scaling process, scaling can be readily achieved, errors 
easily detected, and any further scaling rapidly completed. 
In addition, a completely general analog computer pro- 
gram can be obtained by symbolically retaining the maxi- 
mum magnitude of the problem variables rather than 
computing their numerical values. The following exam- 
ple illustrates the use of this seven-step procedure for 
achieving unity scaling. 


An Example. Consider a hydraulic transient problem in- 
volving two tanks as illustrated in Fig. 6. 

Assume that the cross-sectional area of tank 1 is A, 
and that of tank 2 is A,. Further assume that the volu- 
metric flow of fluid through the valves is linearly related 
to the difference in liquid level across the valves. Given 
initial values for h; and hy, it is desired to determine the 
variation in these levels with time. The equations which 
describe this system are given as follows: 


For tank 1 
A, (dh,/dt) = 4, — 45 (6) 
For tank 2 
A, (dh,/dt) =9,+93— 4% (7) 
For resistance 1 
Ig = (hy — h)/R, (8) 
For resistance 2 
4 = hy/Ry (9) 
Or 
(10) 
And 


(1) 
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Fig. 7—This analog circuit represents the variables in the 
two-tank example problem. 


with initial conditions of h, (0) and h, (0) respectively. 

Assuming that q, and q» are positive constants, a pre- 
liminary analog diagram for solving Equations (10) and 
(11) for h, and hy versus time can be drawn. Prior to 
specification of numerical values, the diagram in Fig. 7 
and the list of variables which are outputs from amplifiers 
and the pot settings in Table 2 are obtained. 

Letting a superscript* indicate the maximum value of 
the output from the amplifiers as indicated in Table 2, 
the process of unity scaling can be symbolically carried 
to completion. The results obtained are given in the final 
column of Table 2. The output from each amplifier be- 
comes the value of the problem variable divided by its 
maximum value and some of the values of the pots are 
changed to compensate in the change of scale from one 
amplifier to the next. The latter changes perhaps become 
more evident if the final analog diagram given in Fig. 8 
is considered. 


In readjusting pot values following the scaling of out- 
puts from amplifiers, four cases arise: 


Case 1. The pot value must be divided by a maximum 
expected value. Examples are pots 1, 2, 4 and 5 wherein 
the maximum value is included so that the inputs to the 
amplifier contain the same scale factor as do the outputs. 


Case 2, The pot value remains unchanged. Examples 
are pots 3 and 6 which remain unchanged because the 
input to the pot contains the same scale factor as is de- 
sired in the output. 


Case 3. The pot value is changed by a ratio of expected 
maximum values. An example is pot 7 which is changed 
by /is*/h,* because the input to the pot contains the fac- 
tor 1/h,* whereas the output must contain the factor 
1/h,*. Pot 8 is another similar example. 


Case 4. Additional pots must be incorporated into the 
circuit, This would arise, for example, had h, — hz been 
desired. In this case inputs of —h,/h,* and hz/h.* would 
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Fig. 8—Amplitude scaling determines the proper setting for 
the amplifiers and potentiometers. 


have been available and an output of (hy — hz) /(/t — ha) *, 
where (h, — h2)* #h,* or h:*, would have been the 
scaled output from the summer. Additional pots, not 
found necessary in drawing a preliminary diagram, would 
have been needed in the input lines to compensate for 
the change in scale factors. 

Choosing the numerical values given in Table 3, a nu- 
merical solution can be obtained on the analog computer 
once h,* and h,* are estimated. These may be estimated 
by calculating the final values of hy and he at equilibrium, 
noting that these values are larger than the initial condi- 
tions, and conservatively choosing maximums which are 
about twice the steady state values to account for the 
possibility of oscillation before the final values are at- 
tained. The numerical values of the variables and the pot 
settings for this numerical example are given in Table 4. 

It should be noted that the analog diagram for the 
foregoing example does not include physical constraints 
which may be imposed on the actual situation. Conse- 
quently, for certain parameter values, the analog solu- 


TABLE 3—Parameters Used for the Example Problenr™ ww" 


SSS 


Physical Parameter Numerical Value 





A, 2 sq. ft. 
A, 45sq. ft. 
% 8 cu. ft./min. 
% 6 cu. ft,/min. 
R, 2 min./sq. ferme 
R, 1 min,/sq, ft. 
h,(0) 15 ft. 
h,(0) 10 ft. 
h, at equilibrium 30 ft. 
h,* = 60 ft. ae 
hy at equilibrium 14 ft. 
hy* ce 30 ft. 


— See 
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TABLE 4—Settings Used for the Example Problem 








Component Numerical Value 
Amp 1 h,/60 
Amp 2 —h,/60 
Amp 3 h,/30 
Amp 4 —h,/30 
Pot 1 0.2500 
Pot 2. (through gain of 0.1) 0.6667 
Pot 3 0.2500 
Pot 4 0.3333 

Pot 5 (through gain of 0.1) 0.5000 
Pot 6 0.3750 
Pot 7 0.1250 
Pot 8 0.2500 





tion will predict negative values of h; and/or hy. In a 
similar vein, values of h; and/or hy may be obtained 
which would represent an overflowing of the tanks, The 
manner in which such constraints can be built into an 
analog program will be discussed in a later article. 


The maximum magnitudes for this example were esti- 
mated from a knowledge of the behavior of the physical 
situation. In other cases, one must resort to other tech- 
niques for estimating maximum and minimum values. 
One method is to obtain these values by a trial com- 
puter run. Another technique is that discussed by Jack- 
son’ which is called the equal-coefficient rule. For differ- 
ential equations of second order or higher, this rule may 
frequently be used to estimate the maximum magnitude 
of the dependent variable and of its derivatives. Given 
an n-th order differential equation of the form 


avy d¥-1y 


dy = 
ay Sy Fava gaa tee ha Ge taht) (12) 


where f(t) is the forcing function, the maximum magni- 
tude of y, 

dy d¥y 

anh? Ds Stine 


are estimated as follows: 


Replace f(t) by its maximum magnitude, denoted by 
A. The dependent variable and its derivatives are multi- 
plied and divided by their maximum magnitudes to yield 





» [ ay 
Ny a ARES 
sees eee = 13. 
ay| ar Fixe +. + dy [y] pi |=4 (13) 
ty 
Having numerical values for A, a, a, ..., ay, the 
maximum magnitudes are estimated by assuming 
ay 1* av-ayy* a]* 
ov[ $e Fees [ See a rane ss 





REPRINTED FROM HYDROCARBON PROCESSING 


and a, [y]* =24. (15) 


Three precautions should be exercised when using this 
rule. 


1. The estimated maximum magnitudes must either 
continually decrease or increase as the derivatives are 
considered in order. If this is not the case, the magni- 
tudes should be estimated in another manner. 


2. The rule assumes zero initial conditions. For a prob- 
lem with non-zero initial conditions, the estimated 
magnitudes should be compared with the initial con- 
ditions and the larger of the two in each case used 
for scaling purposes. 


3. The rule assumes that the problem has a stable so- 
lution. If the programer is doubtful of his problem’s 
stability, the stability is easily examined by using the 
Routh criterion. 


As an example of the equal-coefficient rule consider 
the hydraulic transient problem and the determination 
of h,*. Combining Equations “” and “”), eliminating the 
dependence of fh; on hz, the following equation is ob- 
tained. 


ah, 1 1 1 dh, 1 we 
ca | ape temn, teasR; | cde ADR, = 
1 1 1 
eee Peen of we 1 
stgelete|s © 0M 
Using the equal coefficient rule and evaluating using the 


numerical values in Table 3: 
(h]* =2 (aR +1 





1 
A,R,A, 


R, + 4,R,) = 60 (17) 

















q2 
dh,|* _ AR Sar (1s) 
| ee 
aR, AR, 
qe a cr 
= ie ease = t475 19) 
A,R,A, A,A,R, A,A,R, ‘ , 





The magnitudes decrease as the order of the derivative 
increases, a check of the initial conditions indicate that 
they are less than the estimates, and the solution is 
known to be stable. Hence the indicated values may be 
used as estimates for the maximum magnitudes. 

In many problems the computer operating time, in 
addition to the output voltage from amplifiers, must be 
scaled. If the physical problem takes hours or days for 
completion, the analog solution must be speeded up both 
for convenience and precision. Likewise if the physical 
problem is completed very rapidly, the analog solution 
must be slowed down, In the next article of this series, 
procedures which enable time scaling to be achieved will 
be presented. 
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PART 3: TIME SCALING—When the 
actual time for a process change is 
very short or very long, time scaling is 
used to make it better suited to 
analog computer speed 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


THERE ARE A LARGE NUMBER of physical problems in 
which the variable of interest (the dependent variable) 
changes either very rapidly or very slowly. Reaction kinet- 
ics problems in which the concentrations of some of the 
components change in fractions of seconds are often 
found in the hydrocarbon processing industry. The com- 
bustion of a hydrocarbon is an example of a problem in 
which it may be very difficult to determine the reaction 
rate or follow concentration with time by using a com- 
puter running at the same speed as the problem (real 
time). For such a problem, voltage changes occurring in 
the computer may be so rapid that some of the compu- 
tational elements may not accurately follow the changes 
because of the limited precision of analog components 
for high frequency voltage fluctuations. 

Problems which require a great deal of time to ex- 
hibit the phenomena of interest are also important to 
practicing engineers. The study of the dynamics of large 
scale process equipment falls in this category. Because of 
the excessive amount of computer time which may be 
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involved and the effect of noise and amplifier drift on 
a solution’s validity, the analog simulation of such pro- 
cesses in real time may not be satisfactory. The ability 
to either speed up or slow down the solution of analog 
simulation is one of the important characteristics of an- 
alog computation. The process of scaling the indepen- 
dent variable is called time scaling. 


Recording The Solution. Quite often it is desirable to 
record the results of an analog computation. The prob- 
lem frequencies must then be adjusted by time scaling 
so that the recording device can follow the results. There 
are in general three types of recording devices commonly 
used to record analog results; cathode-ray oscilloscopes, 
galvanometer-type recorders, and servo-motor driven re- 
corders, Cathode-ray oscilloscopes with ac and de ampli- 
fiers can be used for the entire frequency range encoun- 
tered in analog computations, while galvanometer-type 
recorders should not normally be used at frequencies 
greater than 60 cycles/sec., and seryo-motor driven re- 
corders, such as X-Y recorders, are restricted to frequen- 
cies less than about 2 cycles/sec. 





C3/Cy 


Ca fC 
-Ref. | 


= 


CofC 














dy/dt {] -y +> 











Fig. 9—The differential equation as depicted by an unscaled 
circuit. 


LEARN ABOUT ANALOG COMPUTERS . . . 





3/1 p* 








C4/C,p? 
Ret | 





dy/dr ] +y 


[> 


C2/CiB 








Fig. 10—The circuit after time scaling is considered. 























Fig. 11—The time-scaled circuit may be rearranged. 


Solution Oscillation Frequency. It is often helpful be- 
fore attempting to program a problem on the analog 
computer to make an estimate of the frequency of oscil- 
lation of the solution so that a decision can be made as 
to whether the problem will need time scaling. The fre- 
quency of oscillation of the solution of simple differential 
equations can be readily estimated. For a first-order differ- 
ential equation the reciprocal of the time constant is a 
rough estimate of the problem frequency. For example 
for: 


C, (dy/d®) + Cyy = Cyx 
T=1/e=C,/C, 


For this type of equation the solution of the problem will 
have progressed to within 63 percent of its final value 
within a period of time equal to the time constant T. 


TABLE 5—An Example Problem 





Component Output Max Value Scaled Value 
Amp 1 a a, a/a, 
Amp 3 6 Un b/b, 
Amp 5 c oy c/e, 

Pot 9 1 = a/c, 

Pot 10 1 = b/c, 
Pot 1 a(0) a, a(0)/a, 
Pot 2 ky = (b,/a, ke 
Pot 3 k, = ky 

Pot 4 0 _ 0 

Pot 5 ky ES (a,/b,)k, 
Pot 6 Ke E (c,/b, )ky 
Pot 7 ky + ky = (ky + ky) 
Pot 8 n n n/c, 
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For a second order differential equation of the form 
C,(d*y/d9*) + C2(dy/de) + Czy = C,x; the natural fre- 
quency of the equation is o, =(C;/C,)*. The natural 
frequency, for a third order differential equation of the 
form C,(d°y/de*) + C2(d?y/de?) + C;(dy/de) + Cyy = 
Csx; is on =(C,/C2)*. There are no simple approxima- 
tions for determining the natural frequency of differen- 
tial equations higher than third order. However, if the 
coefficients of the low-order terms of an equation are 
quite different from the coefficients of the high-order 
terms, a change of time scale will probably be required. 


Differential Equations. When it is found that a change 
in time scale is desirable for a differential equation, a sub- 
stitution must be made. If t is the real time variable, and 
+ is the computer time variable, the time for solution of 
the equation may be changed by substituting +/f for t. 
If the solution is to be slowed down, 8 > 1; and if the 
solution is to be speeded up, 8 < 1. As an example we 
will change the time scale of the following equation by 
the factor B. 


C,(d?y/dt®) + C, (dy/dt) + Csy=C, (20) 


We substitute 7/8 for t. 
C, [d/d(r/B)} [dy/d(7/B)] + C, (dy/d(r/B)] + Cy = Cy 
C,B? (d*y/d7?) + CyB (dy/dt) +Cyy=C, (21) 

The undamped natural frequency for the unscaled 
equation is (C;/C,)% and for the time scaled equation 
it is (C/ (CB?) )* or (1/8) (Cs/C1) *. Thus, the rate at 
which the solution is obtained has been changed by the 
factor 1/8. It should be noted that 8 appears only in 
terms containing time derivatives and in the forcing func- 
tion if it is time dependent. 

It is of interest here to see how time scaling of the 
last relation modifies the computer circuit. A circuit 
for the unscaled Equation (20) is shown in Fig. 9. The 
time-scaled Equation (21) is shown in Fig. 10. The cir- 
cuit may then be rearranged to give Fig. 11. 

The form of Fig. 11 is particularly significant. It should 
be noted that the circuit is identical to the unscaled cir- 
cuit of Fig. 9 with the exception that the input to each 
integrator is scaled by the factor 1/8, the time scaling 
factor. Thus, the process of time scaling is much simpler 
than magnitude scaling (See Part 2) and in general in- 


TABLE 6—How Circuit Components Are Scaled 





Physical Parameter Numerical Value Scaled Value 


Amp 1 [0.014] [0.014] 
Amp 2 [0.015] (0.015) 
Amp 5 [0.01c] [0.01¢] 
Pot 1 1.00 1.00 
Pot 2 50 0.50 
Pot 3 100 1,00 
Pot 4 0 0 

Pot 5 100 1,00 
Pot 6 5 0.75 
Pot 7 200 2.00 
Pot 8 1.00 1.00 











volves changing 
factor 1/B. 


ator by the 


Methods of Time Scaling. There are two general meth- 
ods of time scaling a problem. One method involves the 
formal time scaling of the differential equations at the 
same time that the equations are magnitude scaled. The 
second method involves magnitude scaling the problem 
and subsequently time scaling the problem by changing 
each of the input pots to integrators by 1/8 after the 
problem has been patched. The second method is usually 
recommended because it permits the magnitude and time 
scaling processes to be done separately and therefore is 
not as confusing as the first method. 

If after magnitude scaling, it is found that pots to 
integrators have settings of less than 0.1 or greater than 
30, time scaling will usually be necessary. In some cases 
after amplitude scaling a situation such as that shown 
in Fig. 12 arises, with the magnitude of other pots to 
integrators about 0.0001 and all other pots in the circuit 
between 0.1 and 30. A time scale factor of B about 
0.0005 is required, for all the integrators except the in- 
tegrator shown in the figure above. The pots which are 
inputs to the summer may be scaled by 1/8 instead of 
the input pot to the integrator. The output from the 
integrator is then time scaled but the output of the sum- 
mer contains a 1/f factor. 


The steps in time scaling can be reduced in many 
cases to: 


Step 1. Changing the input of each integrator by 1/8. 
> 1 slows the solution down, 8 <1 speeds the prob- 
lem up. 


Step 2. Choose 8 such that pot settings are between 
0.1 and 30. Settings close to 1 are the most desired. 


Step 3. Change by a factor of 0.5 or 2 and replot 
the solution. If the curves are identical the problem is 
properly time scaled. 


Rate of Integration. It is of interest to see how inte- 
gration is slowed down or speeded up by changing the 
magnitude of an input signal to an integrator. A simple 
integration circuit is given in Fig. 13. 

The output voltage is 


=—(1/RC)5S,¢e,dt+ E, 





where Ey is the value of eo at t= 0. The group 1/RC 
has the units of reciprocal time and along with the 
voltage ei controls the rate at which €. changes. In many 
computers R = 10° ohms and C = 10° farads so that 
RC normally equals 1 second. There is also usually a 
means for changing the capacitor C by a factor of 10. 
The resistance R can of course be varied by inserting a 
pot in the input of an integrator. If RC is one second 
then the rate at which the output voltage of the inte- 
grator changes is just ¢;. Now if RC is greater than one 
second corresponding to 8 > | the rate of change of é 
will be decreased, and if RC is less than one second cor- 
responding to 8 < 1 the rate of change of é, will increase. 


An Example Problem. To illustrate the rules for time 
scaling we consider a constant volume, isothermal, gase- 
ous reaction represented by the following relation. 


16 














ei { & 


























Fig. 14—A circuit to solve the simultaneous equations. 


ky ks 
A2BeC (22) 
ky ky 


We will let a, 6 and c be the moles of A, B and C re- 
spectively, at any instant and ka, ke, ks and ky will be 
the rate constants expressed in reciprocal seconds for the 
various reactions shown in Equation (22). 

The rate equations for the change of A, B and C with 
time are: 


da/dt =—k,a+k,b (23) 
db/dt =k,a—k,b—kyb +k, (24) 
de/dt =k,b—kye (25) 


Another relation expressing the fact that the total num- 
ber of moles n remains constant may be written although 
the relation does not introduce any information not con- 
tained in Equations (23), (24) and (25). 
atb+e (26) 

If it is desired to follow the concentration changes of 
a, b and c with time, Equations (23), (24) and (26) 


n 
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Fig. 15—Summation can be performed in the integrators. 
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Fig. 16—After amplitude and time scaling is applied. 


uF 
Ref —o—[>— Ramp 


Fig. 17—The time signal for the plotter is also scaled. 


must be programed on the computer. The physical param- 
eters chosen for illustration are given below. 








ky = 100 sec"? 
k, = 50 sect 
ky = 150 sect 
75 sect 
100 
0 
(0) =0 


We may now draw a preliminary computer circuit to 
solve these simultaneous equations as shown in Fig. 14. 

This circuit may be simplified if we are not interested 
in obtaining da/dt or db/dt by performing the summa- 
tion in the integrator. The revised circuit is shown in Fig. 
15. 

We now must scale Equations (23), (24) and (26). 
We will first magnitude scale and then time scale if neces- 
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sary. It should be noted that a, b and ¢ each must be 
magnitude scaled. The reactions are constrained and the 
initial conditions are such that the concentrations of 
A, B or C at any time may not exceed 100, the initial 
concentration of A. 

The numerical values of each of the physical param- 
eters may now be determined, assuming a, = b; = ¢, = 
100. These are shown in the second column of Table 6. 

The large gains on the integrators indicate that the sys- 
tem is a naturally fast one and must be time-scaled to 
slow it down. Since it is desirable to have pot settings 
beween about 30 and 0.1 we choose the time scale factor 
B to be 100. 

Only pots that are on inputs (but not initial condition 
pots) to integrators must be scaled. Pots 2, 3, 5, 6 and 7 
are adjusted by dividing the numerical value found from 
magnitude scaling by 8 = 100. The magnitude and time 
scaled values are listed in the third column of Table 6. 
The final computer circuit (amplitude and time scaled) is 
given in Fig. 16. 

If it is desirable to plot the concentration of A, B, or C 
as a function of time on an X-Y plotter, a time signal 
for the plotter may be supplied by the analog computer. 
It is important to note that the analog circuit for the time 
signal (a ramp increase in voltage) should be scaled just 
as the problem. Suppose the solution is to be recorded in 
F seconds of computer operating time. The scaled circuit 
for the ramp signal is shown in Fig. 17. 

The ramp units are equal to 7/F, where 7 is the com- 
puter operating time, and thus will reach a maximum 
value equal to that of the reference voltage after F sec- 
onds of operating time. 

In order to obtain the relationship between actual 
problem time t and the length of computer chart travel, 
two additional factors must be known: the time scale fac- 
tor, 8 =7/t; and the recorder chart setting, G, volts per 
inch. The relationship is developed in the following man- 
ner: 

Ramp = (Ref. Volts) +/F, volts 
Since + = Bt, 

Volts Developed = (Ref. Volts) Bt/F 

Volts Developed per Unit Actual Problem Time 

= (Ref. Volts) B/F 
Using a recorder chart setting of G volts per inch, 

Inches of Chart Swept by Ramp per Unit Actual Prob- 

lem Time = (Ref. Volts) B/(FG) 


As an example, if the solution is to be recorded for 
10 seconds of computer operating time using a 10-volt 
machine, a recorder set at 1 volt/inch, and 8 = 10, gives 
10 inches of chart = 1 unit actual problem time. 

In many problems it is desirable to check to see if the 
chosen time scale is such that the solution is recorded 
accurately. This may be done by plotting the solution for 
the chosen time scale and then replot the solution for a 
different value of 8. If the two solutions do not overlap 
the chosen time scale is a poor one and some of the solu- 
tion information is being lost. New time scales should be 
chosen until two different values of f give solutions 
which coincide. 


Indexing Terms: Analogs-9, Circuits-10, Computations-4, Computers-9, De- 
scriptions-8, Electricity-10, Engineering-4, Programing-10, Simulation-4, 
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PART 4: PROCESS SIMULATION-There 
is a usual sequence of steps which are 
required to complete an analog 
simulation 


Theodore G. Smith and Theodore W. Cadman 
University of Maryland, College Park, Md. 


Now we SHALL piscuss the steps before patching, 
standard modes of analog operation, suggestions to fa- 
cilitate accurate patching, and the operation of the ana- 
log together with the interpretation of results, The aspects 
of analog simulation are considered first in general terms 
and then as they apply to the simulation of the vapor- 
phase dehydrogenation of benzene in a homogenous, iso- 
thermal, flow reactor. 


Steps Before Patching. In order to adequately prepare 
a problem for analog simulation, four principal steps 
should be carried to as near completion as possible be- 
fore patching is undertaken. These steps are: 


1. The development of the mathematical model for 
the system to be simulated 

2. The preparation of an analog diagram exclusive of 
scaling 

3. Amplitude scaling of the analog circuit, and 

4. Time scaling of the analog circuit. 


Although the first step is an obvious one, it is perhaps 
the most crucial for all subsequent aspects of the simu- 
lation largely depend upon the results of this first step. 
Apart from the difficulties, familiar to all engineers, en- 
countered in mathematically describing a physical sys- 
tem, several mathematically equivalent forms can 
frequently be obtained for a given system. Since the 
relative ease of incorporating the various forms as well 
as the accuracy of the analog simulation may differ, it 
is wise to capitalize on the advantages of the analog and 
to minimize potential analog trouble spots while develop- 
ing the model. This is perhaps more readily done by 
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considering Step 1 and Step 2 simultaneously even if 
the latter step is only partially completed. Further con- 
sideration will be given to analog trouble spots in the 
next article of this series. 

As an example of how a knowledge of analog char- 
acteristics can simplify the analog circuit and inhance 
e—** cos wt 
where @ and w are constant. In this form the generation 


accuracy, consider the generation of y= 





of y would require an integrator to produce t, an expon- 
ential circuit to generate e—*', a cosine generator, and a 
multiplier to complete the simulation. Three special ana- 
log components, generally of very limited availability, are 
required and nonlinear components are connected in 
series thus reducing analog accuracy. Recalling that the 
analog permits accurate and continuous integration with 
respect to one independent variable at almost any speed, 
consider the differential equation 


d*y/dt® + 2a dy/dt + a? + w* = 0, 


with initial conditions of 1 and —a for y and dy/dt re- 
spectively. Its solution y = e—*' cos wt is seen to result in 
a simpler and more’ accurate analog circuit. Two inte- 
grators and an inverter are the only major components 
required. 

Frequently of equal or greater concern than ease and 
accuracy in simulation is flexibility in the analog circuit 
because the analog can be an excellent tool for rapidly 
examining the effect of changes in process parameters 
and characteristics. Flexibility arises from retaining a 
one-to-one correspondence between the physical vari- 
ables and the analog variables. Thus although an elec- 
tronic analog computer circuit is not an electrical ana- 
log of the physical system in the conventional sense 
(since all physical variables are represented by analog 
voltages), it should be considered as more than a group 
of components designed to perform specified mathematical 
operations. An analog circuit should represent flexible 
simulation of the physical system. Such a circuit, which is 
called an analog by many, is obtainable through the one- 
to-one correspondence of physical and computer variables. 
The degree of flexibility to be built in an analog circuit 
is, of course, dependent upon the modifications which may 
be made and this should be seriously considered during 
problem definition. 
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As a simple example of flexibility, consider the con- 
secutive reaction steps described by 


Ps 15, pg’ ¢ 
dA/dt = —hA (27) 
dB/dt = kA — ky B (28) 
dC/dt = ky B (29) 


A preliminary analog circuit based on these equations is 
given in Fig. 18. An alternate, but equivalent, solution 
can be obtained by eliminating the dependent composi- 
tion B to yield: 


dA/dt = —hA (30) 
dC /dt® + k, (dC/dt) — ky k, A=0 (31) 


The preliminary analog diagram for Equations 30 and 
31 is shown in Fig. 19. Consider now the relative ease 
of changing k, B to ky B in Figs. 18 and 19. In the first 
case only a squaring card needs be added after the in- 
tegrator for B. In the latter case, however, it is necessary 
to revert to the original equations and repeat the elimina- 
tion of B. This yields 


a eS 
the incorporation of which would require major modifica- 
tion of Fig. 19. Similar difficulties would be encountered 
for other model changes. Thus, although the circuits in 
both figures represent the same physical system, Fig. 18, 
in which the one-to-one correspondence is retained, is 
significantly more flexible. 

Once the model has been completed, with due consid- 
eration to the simplicity, accuracy, and flexibility of the 
analog circuit, the details of the preliminary analog cir- 
cuit can be completed and the scaling undertaken. In 
completing the scaled diagram, it is suggested that: 
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Fig. 18—Flexibility must be considered 
when setting up an analog computer 
simulation circuit. 


Fig. 19—Although this cir- 
cuit is simpler, it is not 
as flexible to changes. 











ki 
A 
kaky 
-dC/dt ¢ 
ko 


a. Each analog variable be labeled as to its physical 
meaning at each step. This will be of aid not only in 
the scaling process but in the remaining steps of a 
given simulation as well. 


b. The process of unity scaling described in Part 2 of 
this series is recommended for completing amplitude 
scaling. If maximum magnitudes are unknown or may 
change from one run of the final simulation to another, 
amplitude scaling may be completed in a symbolic 
manner so that modifications may be readily made. 


c. The process of time scaling as described in Part 3 of 
this series is a systematic way of completing this step. 
It is to be noted that time scaling does not affect am- 
plitude scaling but that changes in amplitude scaling 
may affect the choice of an appropriate time scale fac- 
tor. Thus time scaling should be completed after am- 
plitude scaling. 
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Modes of Analog Operation. Before discussing the re- 
mainder of the steps generally involved in obtaining an 
analog simulation, it is convenient to describe the modes 
of analog operation. Most analog computers have the 
following modes of operation which can be chosen by the 
computer operator: 


Pot Set Mode. When the computer is in the pot set 
mode, the values of the pots may be set. A reference 
voltage is applied to the input to the pot (this may be 
done automatically) and the value of the pot is set by 
monitoring the output voltage of the pot using a volt- 
meter or a digital voltmeter. The inputs to the amplifiers 
are disconnected and grounded and the amplifiers are 
effectively at zero gain. The analog computer is generally 
placed in the pot set mode when a simulation is not being 
carried out. 


Reset Mode. In the reset mode, the initial conditions 
are imposed on the integrators. The inputs to the inte- 
grators are disconnected but the analog circuit is other- 
wise operable. 


Operate Mode. In the operate mode, the inputs to all 
amplifiers are connected and the initial condition volt- 
ages to the integrators are disconnected. Problem solu- 
tion is affected according to the mathematical operations 
specified by the analog circuit. 


Repetitive Operation Mode. The purpose of the repeti- 
tive operation mode is to permit analog solutions to be 
obtained very rapidly. (8 is typically changed by a fac- 
tor of 0.02 to 0.002 corresponding to a speed up of 50 to 
500.) It enables the display of an analog solution on an 
oscilliscope or the performance of calculations in which 
a rapid speed of solution is desired. In this mode, an 
internal change in the time scale is automatically made 
so that the simulation may be speeded up. The computer 
automatically cycles between the reset and operate 
modes, being in each mode for a specified amount of 
time. The solution obtained during each operating cycle 
can be displayed on an oscilliscope. 


Suggestions for Accurate Patching. Patching is the 
term used to describe the process of connecting together 
the various analog components making up an analog cir- 
cuit by means of external leads. Since the details vary 
from one type of analog to another and are adequately 
described by the manufacturer of each, they are not dis- 
cussed here. The patching process is generally completed 
while the computer is in the pot set mode or while the 
removable patch panel is detached from the analog. Once 
the patching is complete and the patch panel in place on 
the computer, the pots may be set. Since the resistance 
of the pots will vary with the loading, it is important to 
set the pots after patching is complete and when the 
computer is in the pot set mode. (In this mode, the 
inputs to the amplifiers are grounded so that loading of 
the pots occurs.) Because of loading, the voltage output 
from the pot should be used for setting the value rather 
than the scale which may be indicated on the pot. 


Although the analog circuit is now ready for use, the 
patching, the pot settings, and the operability of the 
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components in use should be independently checked to 
eliminate solution errors which may arise from these 
sources. These checks are characteristically made by per- 
forming: 


e@ A static check, and 
@ A dynamic check. 


In completing a static check, arbitrary initial condi- 
tions are placed on all integrators. The outputs which 
would be obtained from all analog components if the 
initial conditions were placed on the integrators but the 
inputs to the integrators were removed are calculated 
from the mathematical model and/or the scaled analog 
diagram. The computer is then placed in the reset mode 
and the outputs obtained from the analog components 
by means of a voltmeter or a digital voltmeter are com- 
pared with the values calculated. In this manner the pot 
values (except for those which are inputs to integrators) 
and the patching may be checked and the static opera- 
bility of the components established. The pots on the 
input lines to integrators can be checked separately. 

The dynamic check corresponds to trial computer runs 
in which the dynamic operability of the analog compon- 
ents is established and the feasibility of the scale factors 
which have been assumed is checked. A computer run is 
carried out by first placing the computer in the reset 
mode to establish the initial conditions on the integrators 
and then in the operate mode, or simply placing the com- 
puter in repetitive operation mode. The computer is left 
in the operate mode for the desired time of the simulation. 


Computer Operation. Once the static and dynamic 
checks have been made and errors corrected, computer 
runs can then be made for the purpose of recording the 
results. If only qualitative results are desired, they may 
be observed on an oscilliscope by placing the computer 
in repetitive operation mode. If a permanent record is 
desired, they may be obtained using an x-y recorder. 
Since the response time of a recorder is quite slow com- 
pared to that of an oscilliscope, the normal speed of 
computer operation must be used when recording, The 
computer is first placed in reset mode to establish initial 
conditions and then in the operate mode during which 
time the results are recorded. 


The accuracy of a problem solution obtained on an 
analog computer is dependent upon several factors. Solu- 
tion accuracy depends upon how faithfully the algebraic 
and differential equations represent the physical system, 
how accurately the physical constants are known, and 
upon how well the problem is programed so that each 
computer component is used in its most accurate range. 
Under the best of circumstances one might expect the 
computer solution to be accurate to at least three sig- 
nificant figures. 

It should be noted that the solution or portions of the 
solution obtained on an analog computer may not be 
obtainable in a physical system. To illustrate: The dif- 
ferential equations describing liquid level in a tank as 
normally solved on an analog computer can yield nega- 
tive heights. This of course is physically impossible. It is, 
therefore, important to have a physical feeling for the 
problem limitations when evaluating computer results. 
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Fig. 20—This preliminary analog diagram retains one-to-one correspondence with the dehydrogenation problem. 


Thus far in describing the steps of a simulation, atten- 
tion has been directed toward the calculation of the 
values of known variables for a specified known system 
with specified known parameter values. This has per- 
mitted the sequential examination of 


© The development of the mathematical model and 
the preliminary analog diagram 

© The development of the scaled analog diagram 

© The achievement of accurate patching, and finally 

© Computer operation leading to the desired results. 


As such, the presentation describes a rather narrow facet 
of analog simulation for frequently the simulation will be 
directed towards: 


a. The calculation of the values of known variables 
for a specified, known system using numerous param- 
eter values 


b. The calculation of the values of known variables 
for numerous modifications of the simulated system 
using numerous parameter values 
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c. The choice modification of a simulated system 

and/or its parameters in order to obtain results possess- 

ing some characteristic, such as the determination of 

the optimum production rate for a chemical reactor 

system. 

The simulations classified under a may be readily 
carried out after the circuit for the first set of parameters 
is prepared for operation by systematically: 


© Setting new parameter values, and 


© Recording the results for each set of parameter 
values. 


If the parameter changes to be made are small, there is 
generally no need to rescale, However, if the changes are 
large, the assumed scaling factors should be periodically 
checked and the circuit rescaled if necessary. In the re- 
scaling process, a clear concise circuit diagram together 
with tables of all analog variables and pot values, which 
presents the results of scaling obtained by symbolically 
retaining maximum values of all amplifier outputs, can 
be of great use. There is generally no need to perform 
additional static and dynamic checks. 
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In completing simulations classified under b, careful 
organization should be begun with the development of 
the initial mathematical model. It is important to build 
as much of the desired flexibility into one analog circuit 
as possible. If this is accomplished, the simulation can 
proceed systematically through 


© Changing the analog circuit, 


© Changing the parameter values over the desired 
range for each circuit, and 


© Recording desired results after each change. 


These last two steps are identical to those employed 
for simulations classified under a. Each change in an 
analog circuit should be viewed as representing the solu- 
tion of a new problem, Despite the frequent temptation 
not to sketch and scale a new diagram or to perform 
static and dynamic checks, the time so spent can repre- 
sent a significant saving over the time which may be 
spent if errors inadvertently creep into the sequence of 
circuit changes. 

Simulations classified under c are the most difficult to 
break down into a series of steps and are frequently sub- 
jected to a more or less random approach because of the 
ease with which circuit parameter changes are made on 
the analog. Such a random approach can result in need- 
less inaccuracies and repetitive examination of the results. 
In order to prevent such occurrences, careful notes of the 
circuit changes made, the parameter values examined, 
and the corresponding characteristics of the results should 
be taken and periodically examined. Scaling should be 
periodically checked when parameter changes are being 
made and each change in patching should be accompa- 
nied by scaling and static and dynamic check, In par- 
ticular, enthusiasm for results should supplement a sys- 
tematic approach rather than replace it. 


Dehydrogenation of Benzene. For purposes of illus- 
trating the steps of an analog simulation, the vapor-phase 
dehydrogenation of benzene in an isothermal, homogen- 
cous, flow reactor has been chosen. The system as de- 
scribed in Smith? using the rate equations reported by 
Hougen and Watson? is used. 

The vapor-phase dehydrogenation of benzene is as- 
sumed to proceed according to 


2 CsHs = CoH + He (33) 
and 
ChHe + CuHio = CisHu + He (34) 


At 1400°F the reaction of Equation 33 gives: 


ry = benzene reacted, Ib. moles/(hr.) (cu. ft.) 


= 6.23 [P, — (PpPx/0.312)] (35) 
and the reaction in Equation 34 gives: 


‘2 = triphenyl produced, Ib. moles/(hr.) (cu. ft.) 
= 3.61 [PsPp — (PrPx/0.480)] (36) 


In Equations 35 and 36, P denotes the partial pressure in 
atmospheres with subscripts B, D, H, and T representing 
benzene, diphenyl, hydrogen, and triphenyl, respectively. 
It is desired to calculate the conversion of an initially 
pure benzene feed at one atmosphere as a function of the 
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Fig. 21—A time ramp is included for recording. 
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Fig. 22—Results for the dehydrogenation simulation. 


Vz/F, ie., the reactor volume divided by the feed flow 
rate. 


Choosing 1 Ib. mole benzene feed as a basis and letting 
x, be the conversion of benzene according to Equation 
33 and x be the conversion of benzene according to 
Equation 34, Smith? shows that 





moles of hydrogen = (3) x1 + x2 (37) 
moles of diphenyl = (4) 1 — x2 (38) 
moles of benzene = 1 — x1 — ¥2 (39) 
moles triphenyl = x (40) 
and that the design equations are 
Va/F = Sdx/n (41) 
and 
Ve/F = Sdx./r2 (42) 


Smith then proceeds to solve Equations 41 and 42 nu- 
merically using the relationships given in Equations 35 
to 40. The method employed is representative of the 
techniques which would be used to simulate such systems 
on a digital computer. 


In order to obtain an analog solution, Equations 41 
and 42 are first expressed in differential form and r, and 
r, are then replaced by the relationships given in Equa- 
tions 35 to 42. This results in: 


dx,/d(Vp/F) = 6.23 
(Q — x1 — x2)* — (x/2 — x2) (1/2 4+ »)/0.312] 
(43) 
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TABLE 7—Unity Scaling, Time Scaling and Static Check 





Amp _ Physical Variable Estimated Maximum 





01 x x* = 05 
02 we x* = 0.1 
03 =x mt = 0.1 
4 Pp=l—-a—m Pz* = 1.0 
05 =Pp a ex) 2 Xe Pp* =02 
06 —Py = —n/2—» Py* = 0.35 
07) ~—Ps' = (1 — x1 —x)° (PsPs)* = 1.0 


08 = —PsPp = 
— (1 — 1 — %) (m1/2 — x2) (PaPp)* = 0.2 
09 PpPy = 
(%1/2 — 2) (x1/2 + x) 
10 PpPy = x2 (x1/2 + x2) 


(PpPx)* = 07 
(PzPz)* = 0.030 





Pot Before Scaling After Amplitude Scaling 








Ol kh = 6.23 (PgPs)* ki/mi* = 12.46 

02. «Ok, 6.23/0.312 (PpPy)* ke/x1* = 2.795 
03 ks = 3.61 (PsPp)* ki/xs* = 7.22 

04 = ky = 3.61/0.480 = (PrPy)* ki/x2* = 2.632 
05 0.5 0.5 */Pp* = 1.25 

06 0.5 0.5 m*/Py* = 0.714 
07 1.0 1.0/P,* 1,000 
08 1.0 x2*/P,* = 0.100 
09 1.0 x1*/P,* = 0.500 
10 1.0 x*/Pp* = 0.500 
11 1.0 x*/Py* = 0.286 


——— a 
Static Check 














Initial Amplifier 
Amp Computer Variable Condition — Output 
01 m/0.5 0.500 0.500 
02 »2/0.1 0.800 0.800 
03 —m/0,1 —0.800 
04 By 0.670 
05 —P,/0.2 —0.225 
06 —Py/0.35 —0.585 
07 —P;° —0.449 
08 —PsPp/0.2 —0.151 
09 P>Py/0.07 0.132 
10 PrPx7/0.035 0.468 
ee eee 

Pot After Amplitude and Time Scaling (8 = 10) 
01 (PaPs)* ki/(x1*8) = 1.246 

02 (PpPx)* ke/(x1*B) = 0.2795 
03 (PaPp)* ks/(x2*8) = 0.722 

04 (PrPu)* ks/(x2*8) = 0.2632 

05 0.5 m*/Pp* = 1.25 

06 0.5 1*/Py* = 0.714 

07 1.0/Pz* = ‘1.000 

08 m*/Pz* = 0.100 

09 m*/Pz* — = 0.500 

10 m*/Pp* = 0.500 

ll x*/Py* = 0.286 


—_).$ 
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dx,/d (Vp/F) = 3.61 
[Q — x1 — x) (x/2 — x2) — (x2) (x1/2 + x4)/0.480] 


44) 
with x = 0 and x: = 0 when V/F = 0. es) 


Retaining a one-to-one correspondence between phys- 
ical and analog variables, the preliminary analog diagram 
given in Fig. 20 is made. Pots 07 to 11 are included be- 
cause they may be required after scaling although not 
required in the preliminary diagram. It is to be noted 
that the independent variable (V/F) is to be computer 
operating time. In Table 7, the physical significance of 
all amplifier outputs and all pot settings is summarized. 

The processes of amplitude and time scaling are then 
sequentially completed as indicated in Table 7. In esti- 
mating the maximum magnitudes, the values of x, and x, 
corresponding to Vz/F = oo were rounded off to the 
nearest 0.1 on the high side ( i.e., the equilibrium values 
were used) ; Pz* was noted to be = 1; P,* and Pp* 
were initially picked using (2) x1* + x.* and (1%) x* 
— x:* but Pp* was changed to 0.2 when a dynamic 
check indicated it to be larger than 0.15; and the maxi- 
mum magnitude for the products were chosen to be 
products of the maximum magnitudes for factors being 
multiplied. 

To complete time scaling, the range of the pot values 
obtained after amplitude scaling were examined. Al- 
though all pot values between 0.1 and 30, the values 
of Pots 01 to 04, namely, those on inputs to integrators, 
were consistently greater than 1 (ranging 2.5 to 12.46) 
indicating that the simulation would be fast. Consequently 
B = 10 was chosen so that a slower simulation could be 
prepared for recording purposes. 

The scaled analog diagram was then drawn. It is the 
same as Fig. 20 except that the amplifier outputs and pot 
values are changed to those quantities given in Table 7. 
Upon completion of the patching, a static check was per- 
formed as indicated in Table 7, A dynamic check was 
subsequently made to check the dynamic operability of 
the components used and the feasibility of the scale fac- 
tors assumed. It was ascertained that the solution was 
completed in about 10 seconds. Following the setting of a 
time ramp for recording as shown in Fig. 21, the results 
using 8 = 10 and 8 = 5 were also recorded and com- 
pared. 

The results shown in Fig. 22 were obtained using a 
scale factor on the recorder of one inch per volt and a 
10-volt analog computer. 

The next part of this series will discuss the generation 
of particular functions using common analog components, 
Consideration will be given to the operational aspects of 
the analog computer in order to assure stable, accurate 
simulation. Of particular concern will be methods for car- 
tying out algebraic solutions using the analog computer. 
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PART 5: CIRCUIT DESIGN-Stability, 
accuracy and flexibility must be 
considered before choosing from 
among several circuit designs 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


Durine THE course of preparing a problem for analog 
solution, several mathematically equivalent methods may 
occur to the programer. At this point, he must make a 
judgment as to which will best suit his particular pur- 
poses. As indicated in the last article of this series, the 
amount of required analog equipment, the accuracy of 
required analog components, and the flexibility of the 
analog circuits can be significant factors in this choosing 
process. Attention can be given to such considerations, 
however, only after analog stability has been assured, 

Below attention is drawn to those situations in which 
analog instability can arise. Particular emphasis is given 
to its occurrence during the solution of algebraic equations 
on the analog. Methods for insuring stability are discussed 
and are extended to indicate the applicability of the ana- 
log in optimization calculations. The article is concluded 
by considering the synthesis of transfer functions on the 
analog. 


Feedback Stability. Although the analog computer per- 
mits the almost instantaneous solution of algebraic equa- 
tions, precautions must be taken in the design of circuits 
to solve algebraic equations, As an example, consider the 
simple equation y = ay + x,¢ > 0, where it is desired 
to solve for y given x. An analog circuit is readily designed 
for the equation in the form given in Fig. 23. 

The solution which might be expected is y=x/(1 — a). 
In actuality, the circuit will become unstable if a > 1 and 
the amplifiers will readily saturate. This results from the 
phenomena known as positive feedback with excessive 
gain, To illustrate the reason for the instability, assume 
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Fig. 23—An easy circuit to solve an algebraic equation al 
though it may be unstable in operation. 











that the output from amplifier 1 is momentarily in error 
by an amount e. Following the progress of the error 
through the feedback loop in an approximate, sequential 
manner, the error ¢ is seen to result in an error of ea in 
the output from amplifier 1 after the first pass. Thus the 
initial error results in a subsequent error of the same sign, 
hence the term positive feedback. Now if a < 1, the error 
will soon die out. However, if a > 1 the error is propa- 
gated and the amplifiers soon overload. Thus if the gain 
ef the feedback loop becomes excessive, a > 1, the circuit 
becomes unstable. 


Positive feedback arises when an even number of am- 
plifiers occur in a loop of an analog circuit, To insure 
stability in such cases, the total gain around the loop 
must be less than unity. 

Negative feedback arises when an odd number of am- 
plifiers occur in a loop of an analog circuit. As indicated 
by W. Brunner,‘ rather restrictive maximum permissible 
gains arise in this case also. These are less than 8, 2.9 
and 2.1 for loops containing 3, 5 and 7 amplifiers respec- 
tively. The presence of integrators in the feedback loop, 
however, markedly enhances the stability, to the point 
where maximum permissible gains are higher than ever 
required. 


To avoid feedback instability it is advisable to: 


© Avoid, where possible, feedback loops which do not 
contain at least one integrator* 

© Avoid summing-amplifier gains in excess of 30° 

© Double check the occurrence of loops containing 
positive feedback (i.e. loops with an even number of am- 
plifiers and integrators to insure that they are required 
by the physical system being examined). 

© Adhere religiously to the maximum permissible feed- 
back gains if loops with only amplifiers are used. 

Such precautions are of particular significance in the 
design of analog circuits for solving algebraic equations. 
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Fig. 24—Solving three equations simultaneously. 














Fig. 25—The upper circuit is a shorthand symbol for the 
more complete circuit shown in the lower portion. 


Algebraic Equations. Algebraic equations appear fre- 
quently in chemical and petroleum engineering as being 
representative of process performance, such as in steady- 
state mass and energy calculations, or as an integral part 
of a dynamic simulation. An example of the latter would 
arise during the simulation of the linear dynamics of mul- 
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ticomponent distillation on an analog. In this case, if 
Raoult’s Law is assumed, it can be shown® that the lin- 
earized equilibrium relationship is 








si ik y 
95,4 = Ky,4 84,4 bay a = 
xe : 
pet 1,8 
where 
¥;, 4 =the deviation from its steady-state value of the mole 


fraction of component i in the vapor leaving tray j 


2), 4 = the steady-state value of the mole fraction of component 


i in the liquid on tray j 

x); =the deviation from its steady-state value of the mole 
fraction of component i in the liquid on tray j 

‘K,, ; = ratio of the vapor pressure of component i evaluated at 
the temperature of tray j to the column pressure 

K’,, ; =the derivative of K;, ; evaluated at the temperature of 
tray j 

N =number of components present. 


To calculate the deviations x;,; given the deviations 
9, i, it is necessary to solve the coupled set of algebraic 
equations. Thus for a four-component system and a given 
tray, equations of the following form arise where the a;,’s 
are constants 





Ya = yy Fy + ae Fa + as (45) 

Va = May Hy + gy Xp + Mpg Xp (46) 

Vy = yy ¥1 gy Hz + gy Xa (47) 
where 

RE 8 

x =—*,—*,— 4, 


One method of analog simulation would be to solve Equa- 
tions 45, 46 and 47 simultaneously to yield 





4 = by I + bia Ia + bis Ms (48) 
= be1 4 t bea Ya + bas Ys (49) 
y= bap 1 + by2 Yo + Bs Is (50) 


and program Equations 48, 49 and 50 directly. Such an 
approach, however, requires considerable computation be- 
fore the analog solution can be obtained and results in a 
relatively inflexible analog diagram since the known co- 
efficients aj;’s are not isolated on individual pots. 

A much more direct approach is to use the analog to 
perform the calculation using the a;;’s directly. This can 
be achieved by rearranging Equations 45, 46 and 47 to 
yield 





= (01 2 2 — 9 Ma)/Oay (51) 
%2= (Yo — 21 ¥1 — M99 %9)/ 420 (52) 
Hy = (Ya — O51 1 — M2 ¥2)/ M5 (53) 


which can be programed in the manner shown in Fig. 24. 


In Fig. 24, the symbol shown on the upper circuit of 
Fig. 25 is used to represent the same as the more com- 
plete components shown on the lower view of Fig. 25, 
where P,, Pz, and Ps are grounded pots. It is readily veri- 
fied that x = — y (P:/Ps) — z (P2/Ps). Thus a factor 
common to all inputs can be handled by the addition of 
the grounded pot P; in the feedback line. 

Also note that Fig. 24 has not been scaled and that all 
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of the aij’s have been assumed to be positive. Negative 
a;;’s are handled by the addition of inverters in the ap- 
propriate feedback lines. 

An examination of Fig. 24 indicates that there are three 
positive feedback loops. 

© The loop consisting of amplifiers 1 and 2 which has 
an over-all gain of (@12 @21) /(@11 @z2) 5 

© The loop consisting of amplifiers 2 and 3 which has 
an over-all gain of (azs @s2) / (@z2 ass), 

© The loop consisting of amplifiers 1 and 3 which has 

an over-all gain of (413 @31) / (11 @ss) - 
If a stable solution is to be obtained each of the three 
gains must be less than unity. For the particular case of 
multicomponent equilibrium, one would generally expect 
ai; > aij; j # i, in which case the stability criteria would 
be automatically satisfied. If, however Equations 45, 46 
and 47 had been solved for x2, xs and x, respectively, the 
resulting circuit would be unstable. 

Although the direct method indicated above can fre- 
quently be used to great advantage, care must be taken 
to insure that maximum permissible feedback gains are 
not exceeded. Particularly troublesome are those cases in 
which the equations can not be rearranged to make 
a;; > aij; j ® i, there is no unique relationship between 
the aj;’s which can be used to the programer’s advantage, 
or the a;;’s are variable. In the latter case, the pots in Fig. 
24 would be replaced by multipliers and the analysis of 
stability could become unduly complex. Under such cir- 
cumstances the programer might be well advised to use a 
method of solution which eliminates the need for a stabil- 
ity analysis. The method of steepest descent is such an 
alternative. 


The method of steepest descent results in the synthesis 
of a set of differential equations, the steady-state solution 
of which represents the solution of the algebraic equa- 
tion(s). The method of synthesis is given below. 

Choosing Equations 45, 46 and 47 for example, errors 
1, & and ey are defined in the following manner. 

8 SI. + ys HH Oye Fe Fas Ms (54) 
Vat yy Ht ae Xp + O29 Xs (55) 
= Ig + 53 *, + O52 Xa + M53 45 (56) 
If the correct values of x:, x2 and x3 are chosen, ¢: = e2 
= e, = 0. Choosing E = *,? + *,? + &,*, E = 0 is seen 
to result if and only if the correct solution is obtained. 

Other forms for E may also be used. For example, 
E = |e,| + |eo| + |es| or E = ex* + ent + 34, B 2 0. 

It is recognized that if an analog circuit is designed for 
which dE/dt < 0 for all t, the steady-state solution would 
result in E = 0 and thus solution. Letting x:, x2 and xs 
be functions of t, the following is obtained by the rules of 
differentiation 


=2 («3 cee «# +6 &) 


&, 








dey de dea oes de. axy 
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From the last form of Equation 57, it is evident that 
dE/dt < 0 is insured if 


oe ae ge 4 428) (08) 
i de de, 
a toa tage) (9) 
dxy _ de de es 
dt -(« a ary ©) 


or in the general case 


Bes Fg OF 
ems Ox; (1) 


For Equations 54, 55 and 56, Equations 58, 59 and 60 
reduce to 


ax,/dt = — (ana + an @ + an 6) (62) 
dxo/dt = — (aw e + am @ + asz 6) (63) 
dxz/dt = — (ay € + ars & + ass €) (64) 


Programing of Equations 62, 63 and 64 as indicated in 
Fig. 26 yields an analog circuit which, after a brief 
period of operation, will produce the desired values 
of x:, x2 and x; without the necessity of performing a 
stability analysis. 


To speed the approach to the steady-state, the speed 
of integration can be increased by using a smaller feed- 
back capacitor. Many analogs permit a 8 = 0.1 or less 
to be selected in this manner. In the application of this 
technique, the calculation of the e;’s is also crucial since 
a null point is being sought. To enhance accuracy, feed- 
back pots may be placed in the feedback line of the am- 
plifiers which calculate the e;’s—Pots P;, P, and Ps in Fig. 
26. By setting these pots at small values, the errors, e;’s, 
generated can be amplified by several hundred so that the 
null point is more clearly defined. If this latter approach 
is taken, momentary overloading of the amplifiers may 
occur prior to the approach to steady-state. It is common 
to limit the values of e;’s between + reference, using 
techniques which are to be described in the next article, 
in order to prevent such overloading. After the analog 
circuit has attained a steady-state for a set of fixed 91, 2 
and ys, variable values of y:, y2 and ys can be handled. 
The values of x,, x2 and x3 obtained from the circuit will 
essentially be the instantaneous solution of the algebraic 
set for the variable y;’s. 

The method of steepest descent as described above for 
linear algebraic equations has several other applications 
and extensions which can be of considerable use to the 
practicing engineer: 
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Fig. 26—The method of steepest descent as a circuit. 


Nonlinear Algebraic and Transcendental Equations. Ap- 
propriate e;’s are defined so that e;’s = 0 corresponds to 
solution. Equation 61 then forms the basis for the syn- 
thesis of differential equations whose steady-state response 
is the solution. If multiple solutions exist, the analog solu- 
tion will proceed in the direction in which dE/dt < 0 for 
all ¢ from the initial conditions to the nearest solution. 


When Number of Equations Exceed the Number of 
Unknowns. For these cases, no change needs to be made 
in the approach. The individual e;’s will, however, reach 
non-zero values at steady-state corresponding to a minimi- 
zation of E. Thus the “best” solution can be obtained ac- 
cording to the relationship chosen between E and the «;’s. 
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The squared error form, E = & e;*, is probably the most 
frequently used. 





Minimization. In light of the application given in the 
preceding method, the extension to the general use of 
minimization becomes a natural one. The approach re- 
mains the same. If E = & e; is to be minimized, Equa- 
tion 61 details the differential equations to be used. If 
E = & ¢; is to be minimized, it may be verified that the 
equations to be used are: 
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If multiple minimums exists, the analog solution will pro- 
ceed in the direction in which dE/dt < 0 for all ¢ from 
the initial conditions to the nearest minimum. 


Maximization. For maximization, the method is re- 
ferred to as the method of steepest ascent. Differential 
equations for which dE/dt > 0 for all t are used. For E = 
% e;?, this is readily achieved deleting the negative sign in 
Equation 61 to yield 


dy _ oe 
di FB |) 
; 


as the required equations. For E = & e;, the equations 
are 


dx; _ y de 
GbE oe x Ox, (67) 
‘ 


For further details of these and other steepest descent 
methods, including optimization with constraints, the 
reader is referred to References 4 and 5. 


Simulation of Transfer Functions. In dealing with 
linear, ordinary, differential equations with constant co- 
efficients, the practicing engineer frequently uses the 
Laplace transform because it permits manipulation and 
solution of the equations using only algebraic techniques. 
(See Reference 7 for details.) Due to the widespread ac- 
ceptance of this transform, particularly in process dynam- 
ics and control, the simulation of transfer functions on the 
analog has received considerable attention. Below several 
relevant aspects of this particular area of analog applica- 
tions are developed. 


Letting s be the Laplace transform, f(s), of a function 
f(t), where f(t) = 0 for t < 0, is given by 


so- [roca 


Based upon this definition, the Laplace transform of stan- 
dard functions and operations can be derived. A short list 
is given in Table 8. 


‘The usefulness of the transform is indicated by consider- 
ing the set of differential equations 


ayjat + y = x,y (0) =0 (68) 
dx/dt +x = z,x(0) =0 (69) 
Upon transforming Equations 68 and 69 they become 
9 () +9 (9) = x (9) (70) 


sx() +x (9) =z) (71) 


which can be algebraically manipulated to yield 
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ay, 


(a) 














Fig. 27—y(s)/x(s) = 1/(s + a). 


IO/ze() = 1/6 +1)" (72) 


where 1/(s + 1)? is the transfer function relating z to y. 
Upon specifying 2(s), y(s) can be determined and y(t) 
recovered using established techniques for inversion. For 
systems involving a number of coupled equations, the use 
of the transform results in considerable saving of time and 
the transfer function becomes a very convenient repre- 
sentation of the system under consideration. 

Consequently, the practicing engineer may frequently 
be concerned with the analog simulation of processes 
which have been characterized by transfer functions 
rather than differential equations. In order to facilitate 
the analog simulation in such cases, techniques for analog 
circuit design directly from the transfer functions have 
been developed. 


Consider for example 
y (s)/x (s) = 1/(8 + a) (73) 


Referring back to Equations 68 and 70, Equation 73 is 
seen to be equivalent to 


avjat+av=x (74) 
the analog solution of which is given in Fig. 27. 
Suppose now that 

J (9)/x (8) = K/\s + a)(s + 4) (75) 


is encountered, By noting that Equations 75 can be equiv- 
alently expressed as 


x (s)/z (9) 


1/(s + a) (76) 
and 


z(9)/x &) 


K/(s + 6) (77) 
the analog circuit is readily designed as given in Fig. 28. 


Adding the acceptable circuit given in Fig 29 for 


TABLE 8—Laplace Transforms 


(0 £9 

A A/s 

t 1/3? 

om 1/(s + a) 
J‘ (O,f0) = 0 sf) 
[soa = 

f @ — T), T constant EF (9) 
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Fig. 28—y(s)/x(s) = K/(s + a) (s + b). 
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Fig. 29—An acceptable circuit for a transfer function. 
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Fig. 30—y(s)/x(s) = 

















= (hs + 1)/(T:s + 1). 
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Fig. 31—y(s)/x(s) = 








= [1 — (1s/2)/[1 + (Ts/2)]. 
2A) 


70) aes soe 1), (78) 


all transfer functions of the form 





= 
(sass Fan) 





ds) 
20> (79) 
can be handled when it is noted from Table 8 that mul- 
tiplication by 1/s is equivalent to another integration. 
Thus with knowledge of the analog circuits for the basic 
transfer functions, the programer can design directly an 
analog program for more complicated transfer functions. 


The basic transfer function 


v(s) _ Tis+1 
x(s) Tas+1 


appears to be more difficult to simulate because the differ- 
ential equation from which this result is 


(80) 
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To (ay/dt) + » = T, (dx/dt) + x (81) 


which implies the need for differentiating the input signal 
x, an operation which is generally unsatisfactory due to 
analog noise. This problem is alleviated by noting that 
Equation 81 is equivalent to 


dz/dt = (x/T:) — (y/T:) (82) 


v= (Ty/T2) x + (Ti/T2) z (83) 


Equations 82 and 83 permit the simulation of Equation 
80 without differentiation as is shown in Fig. 30. 

An approximation to y(s) /x(s) = 71s + 1 is obtained 
by making the value of T/T very large. 

For a more complete listing of analog circuits for stan- 
dard transfer functions than can be given here, the reader 
is referred to References 4 and 8, Note that neither the 
circuits presented here nor in the references are scaled 
and that numerous acceptable circuits may exist for a 
given transfer function. 

The material presented above briefly describes the syn- 
thesis of transfer functions consisting of ratios of poly- 
nomials in s, One of the most important transfer functions 
in chemical and petroleum engineering, however, is 


I(/x() = e™ (84) 


the time delay, which is equivalent to the physical situa- 
tion in which y(t) = x(t — T). This transfer function 
can not be simulated exactly on the analog using only 
standard analog components. In order to approximate the 
time delay, numerous circuits based on expansions of 
have been developed.* *: * The most frequently used is the 
first-order Padé approximation 


2) _'1= Bs? 
x@) 1+ TH/2 (85) 


the circuit for which is given in Fig. 31. 


Other approximations used for greater accuracy are the 
second and fourth Padé approximations and the Stubbs- 
Single approximation. In Reference 8, typical responses 
to step and sinusoidal inputs are given for the various ap- 
proximations. As indicated there, the time delay of a step 
is particularly difficult to approximate. Analog simulations 
using these approximations do, however, represent marked 
improvements over simulations completed by neglecting 
the presence of time delays. 
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PART 6: DIODE SWITCHING-A key 
element in the design of circuits for 
non-linear problems, the diode acts 
as a switch in an operation analogous 
to a check valve 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


PROBLEMS WHICH ARE DISCONTINUOUS in nature are 
often solved on analog computers. Problems with dis- 
continuities can arise in several ways: from a discon- 
tinuous physical system or from the piecewise linear 
approximation of a nonlinear differential equation with 
a set of linear differential equations, each of which is 
applicable over a fixed interval. In the hydrocarbon pro- 
cessing industry it is not unusual to find systems which 
have time delays, hysteresis or backlash, coulomb fric- 
tion, limitations on physical quantities, or are forced 
with discontinuous forcing functions. Each of these non- 
linearities as well as others may be simulated on the 
analog computer. 
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The Diode. One of the key elements used in the circuit 
design of nonlinear elements is the diode. A diode per- 
mits current to flow only in one direction and its op- 
eration is analogous to the operation of a check valve. 
‘There are in general three types of diodes used in analog 
computers. These are the vacuum tube diode, usually 
called a “hot diode,” and two “cold diode” solid state 
devices, the germanium diode and the silicon-junction 
diode often called a Zener diode. The symbols for the 
“hot” and “cold” types are shown in Fig. 32. 


+ HS, 
i P Cc 


Fig. 32—Vacuum tubes (left symbol) and solid state devices 
can be used as switching diodes. 


The latest solid state analog computers use only “cold 
type” diodes which have forward resistances on the 
order of one hundred ohms. If the anode or plate P is 
positive with respect to the cathode C, electrons will flow 
from the cathode to the plate and the direction of cur- 
rent flow will be from the plate to the cathode. If the 
plate is made negative with respect to the cathode, elec- 
trons are repelled, and, as a consequence, no current 
flows. When a diode is conducting, its resistance is a 
nonlinear function of voltage and current. 






Some complex non-linear circuits are possible by combining 
several simpler circuits 


Diodes in Switching. Diodes are used extensively as TABLE 9—Non-linear Analog Computer Diode Circuits 
switching devices in digital computers, but are rarely ; 
used strictly as switches in analog computers. They are _agso.uTe VALUE [eo] 
often used with operational amplifiers and relays to form 


switching circuits. Diodes switch currents rather than iu le = 

voltages and the state of the switch and the magnitude 6 8 
of the current is dependent upon the voltage drop across ' 

the diode. 
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Fig. 33—Minimum selection where x = —min (e,, ey ey - . - 
ey, @,). 
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Diodes are often used in the design of selection and 
comparator circuits. Fig. 33 shows a circuit that may 
be used to obtain the minimum of a set of input volt- =n emt|o 
ages, while Fig. 34 shows a circuit that may be used DEAD SPACE SIMULATION ey 


to obtain the maximum of a set of input voltages. In ‘ 
r% ae 
+0 Ss 


either case, 
on R 
“LR+R, | * 


These circuits require that the source and diode for- 
ward resistance be negligibly small compared with resis- _yverepesig oR BACKLASH 3 
tances R, and Rj. The voltage ¢, may be either positive +2m_lal peo 
or negative. If Ry is equal to Ry and é, is twice the posi- em* [9] Ic 
tive reference voltage, then e, will be the positive refe 
ence voltage, and the circuit in Fig. 33 will find the min’ 1 
mum of the set of input voltages. When e, is positive in : Tal | 
the circuit in Fig. 33 the circuit is equivalent to the AND emt lo] 1 . 
circuit used in digital computers. If e, is negative in the J : 
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circuit in Fig. 34, the circuit is equivalent to the OR 
circuit used in digital computers. 

Diodes may also be used to construct comparator cir- 
cuits. Fig. 35 shows a simple comparator circuit using 
an operational amplifier and two diodes. 

The two ungrounded potentiometers used in the cir- 
cuit have their wipers set at ratios of a and a, of the 
voltages -+ ¢. The potentiometers are used to obtain bias 
for the diodes. The output of the circuit is positive and 
of a value a.¢ when x is less than e2 and is negative and 
of a value ae when x is greater than es. Both x and e: 
may be variables in a problem. 





Relays are often used in conjunction with diodes to 
form analog computer switching circuits. Two types of 
relays are commonly used. They are the differential relay 
of the double throw type, having one or more poles and 
the sensitive-plate relay. Plate-type relays can usually 
be energized by the output of an operational amplifier. 
The circuit shown in Fig. 36 can be used to produce 
the function: 

a ifx-ySo 
t=) ifx-y>o 





The lower diode in the circuit of Fig. 36 limits the 
output of the operational amplifier to —eo volts and pre- 
vents the amplifier from overloading. This circuit is a 
variation of the comparison circuit shown earlier in 
Fig. 35. 
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Fig. 35—A simple comparator circuit using diodes. 
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Fig. 36—A relay circuit for switching. 
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Fig. 37—Complex non-linear circuits can be combined from 
simpler circuits. 


A number of simple circuits have been developed to 
simulate the discontinuous phenomena that arises in 
many physical problems. Table 9 contains some of the 
most often used diode circuits. Each of these circuits 
usually has an equivalent relay circuit replacing the 
diodes, When the diode conducts the relay switch is 
closed. Otherwise, the switch is open. 

Several other diode circuits may be found in the lit- 
erature." 57° 


More Complex Circuits. It is sometimes possible to con- 
struct more complex nonlinear circuits by combining rel- 
atively simple nonlinear circuits in series or parallel. 
Two examples of such circuits are shown in Fig. 37. 


Scaling Nonlinear Circuits. The scaling of nonlinear 
components is accomplished in the same manner as linear 
components. They are scaled by changing the gain of 
integrator inputs. Nonlinear components are affected by 
amplitude scaling of their inputs and it is not unusual 
to have to scale both the input and output. Scaling can 
sometimes be done simply by patching the nonlinear ele- 









Fig. 38—The pressure in the first tank is to be controlled. 
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Fig. 40—Block diagram of the control valve. 








ment on the computer and then scale it while forcing 
the element with a variable input such as a ramp 
function. 


A Process Control Example. As an example of the use 
of a diode circuit in an analog computer problem we 
consider the system of Fig. 38 in which it is desired to 
control the pressure in a tank by employing a pneu- 
matic control valve to throttle the flow. For this example 
the pressure P, in the first tank is to be controlled. The 
symbols are defined as follow 





C = Capacitance, cf./psi 
R = Resistance, psi/cfm 

















Pressure, psi 
Q = Flow rate, cfm 


In this system there are several nonlinear elements. 
Each resistance must usually be considered to be non- 
linear and the pneumatic control valve is nonlinear. For 
this example we will linearize each resistance for small 
departures from the steady state and will consider only 
the control valve as nonlinear. 

A pneumatic valve is usually designed so that it op- 
erates over a pressure range from 3 to 15 psi. If the 
valve is normally closed, it will begin to open when the 
signal from the controller is 3 psi and will be fully open 
when the signal reaches 15 psi. A normally open valve 
operates in the reverse fashion. The simulation of the 
control valve requires the use of a nonlinear element 
(a limiter) because the controller output signal range 
is normally 0 to 20 psi while valve movement is re- 
stricted to the 3 to 15 psi range. 

The outflow Q, from the control valve is a function 
of the upstream pressure P;, the downstream pressure P., 
and the valve stem position X, (or the degree to which 
the valve is open). For this problem we assume that the 
valve is very fast so that valve stem position is linearly 
related (no dynamic element) to the pressure P, supplied 
to the valve bonnet. The relationship between the valve 
elements is shown in Fig. 39. 


The relationship between outflow, pressure drop and 
valve stem position is given in the following linearized 
equation: 


Q(t) = SS x, (y) + 3 Py) + SBP, (y) (66) 





Xy(t) 








Fig. 41—Block diagram of the pressure control system. 
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Fig. 42—An analog diagram for the example problem. 


this in turn may be written as 


Q(t) = oo x, (t) + 20s 4 P(t) (87) 


In the linearized form the partial derivatives are con- 
stants evaluated at the normal operating conditions. The 
block diagram of Fig. 40 illustrates how P,, Pz and X» 
affect outflow Q». 

The relationships between pressure drop, flow rate, re- 
sistance and tank capacitance may be developed. In this 
development it is assumed the tanks are pure capacitors 
(no resistance to flow) and connecting lines are pure 
resistances. 


For Tank 1 
Ci (dP/dt) = Q2 — Q& (88) 
For Tank 2 
C, (dP;/dt) = Qs — Q& (89) 
For Resistance 2 
P, — Py = Qa Re (90) 
For Resistance 3 
Pp— P= QR (91) 





A proportional plus integral controller is used to con- 
trol the process 


P,fe = K.(1 + (1/Td S a (92) 
where P, = controller output pressure 
¢ = error (set point minus feedback signal) 
iy ntegral time 
K, = controller gain 





‘The assumption that the valve is very fast leads to a 
linear relationship between controller output pressure Po 
and valve stem position X. 


Rice KPa (93) 
It should be noted that Equation 93 as it now stands 
indicates that there is no restriction on the extent of 
valve stem movement. The relationship must be modified 
to conform to the physical relation. We have the follow- 
ing restriction on Equation 93. 


i, =n, 3<P,<15 
X,=0 P< 3 
X,=15K, P,> 15 (94) 


The restrictions shown in Equation 94 require a lim- 
iter circuit to be placed in the pressure line before the 
valve bonnet. 
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WITHOUT LIMITER 


WITH LIMITER 








Fig. 43—Valve stem position, —X,. 


WITHOUT LIMITER 


WITH LIMITER 





Fig. 44—Tank pressure, P,. 


In process control problems it is often helpful to draw 
a block diagram of the system. This type of diagram 
shows how signals flow through the system and how the 
various elements are affected by the signals. Fig. 41 is a 
block diagram for the system under discussion. 


An analog circuit may be constructed using the block 
diagram and the limitations on valve stem movement. 
The circuit is given in Fig. 42. 

The normal value of each of the constants in the 
problem is given in Table 10. 

Typical pot settings, determined after magnitude and 
time scaling, are given in Table 11. 

It is now of interest to determine the effect of the 
limiter upon the valve displacement and upon the pres- 
sure P, in the first tank. The disturbance is a step change 
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TABLE 10—Normal Values of Constants 





Function Value 
Sp/10 2.56 psi 
K, 0.7 psi/in. 
Dy 36 min. 
P, 0.036 psi 
Ks 0.044 in./psi 
2O/ex, 10 cfm/in. 
@Q/aaP 0.083 cfm/psi 
R, 45 psi/cim 
Ry 590 psi/cfm 
Papo 3 psi 
Cc, 0.007 cf/psi 
Ge 0.0001 cf/psi 
TABLE 11—Typical Pot Settings 
Pot Function Value 
0 Sp/10 0.023 
1 K, 0.683 
2 V/T; 3,82 
5 Ky 0.44 
6 20/2X, 1.00 
8 2Q/2AP 0.083 
10 1/C, 0.146 
12 1/Ry 2.22 
13 1/C, 0.927 
15 P, (0) 0.038 
20 P, (0) 0.305 
21 P, (0) 0.305 
22 P,/10 0.848 
25 V/R, 0.170 
26 P,/10 0.305 
42 Limiter 0.072 
43 Limiter 0.326 


in pressure caused by a setpoint change. Fig. 43 shows 
how the output from amplifier 2 (in this case —Xy) 
changes as a function of time with the limiter circuit 
and without the limiter circuit. Fig. 44 shows how the 
output of amplifier 10 (in this case P,) changes with 
time with the limiter circuit and without the limiter cir- 
cuit. Because of the nature of the control system (closed 
loop feedback control), the effect of the limiter shown 
in Fig. 43 is not just to truncate the unlimited sinusoidal 
decay curve of Xy. The curve in which the limiter is used 
has been significantly changed in both amplitude and 
frequency and this in turn has changed the pressure 
curve shown in Fig. 44. For this particular case the lim- 
iter has helped to reduce the wild pressure fluctuations 
in the first tank due to a set point change. 


This example has illustrated the importance of de- 
signing the analog circuit to incorporate physical limita- 
tions of the system. The example was simplified to illus- 
trate just the effect of the limiter. A more realistic system 
might take into account nonlinearities in resistances, valve 
dynamics and hysterisis in valve motion. Each of these 
may readily be incorporated in an analog simulation. 
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PART 7: FUNCTION GENERATION— 
Usually the best way to generate 
mathematical functions in an analog 
computer is to simulate differential 
equations. 


Theodore G. Smith and Theodore W. Cadman 
University of Maryland, College Park, Md. 


Ir 1s posstme to generate a function either as the solu- 
tion of a differential equation or by employing a memory 
device which stores the functional relationships. ‘The so- 
lution of differential equations is the strong point of ana- 
log computers, while storage is one of their weakest attri- 
butes. Whenever possible one should try to use the 
differential equation method of function generation. 


Differential Equations. We have referred to differential 
equation function generation in a previous article in this 
series. If we have a function of the form y = f(x), we 
may generate the function as a solution of a n‘* order 
homogeneous differential equation if the function and 
each of its first n derivatives are continuous and if each 
of the derivatives can be calculated. The procedure to 
follow in forming a function by this method is to differ- 
entiate the function a number of times until a homogen- 
eous equation may be formed. The solution of the homo- 
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geneous equation with the proper boundary conditions 
will be the desired function. 

As an example, suppose we required the function y = 
Ke in a portion of a problem. To generate this function 
we differentiate once to obtain 

dy/dt = — aKe™ (95) 
and combine this derivative with the original function to 
obtain a homogeneous differential equation 

(dy/dt) + ay = 0 (96) 
The initial conditions which give the desired solution are 
at t = 0,» = K, and the analog circuit is shown in 
Fig. 45. 


-ref 


y= Ke-at 


Fig. 45—A homogenous differential equation is obtained by 
combining an equation with its derivative. 


If we wish to generate the previous function in a var- 
jable other than machine time we employ a similar 
method. For example, consider the generation of the ex- 
ponential 


y= Ke™ (97) 

whose first derivative is 
dy/dt = — Kae (dx/dt) (98) 

or using symbols gives 
j= — 9% (99) 
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Fig. 46—This differential equation generator requires a deriv- 
ative input. 


The solution of this equation may be obtained from the 
circuit shown in Fig. 46, where yo = Ke~,. The genera- 
tion of this function requires the derivative input x. 


Sine and cosine functions of dependent variables are 
used quite frequently and may be generated as solutions 
of differential equations or by a special device called a 
resolver. If we wish to generate functions of the form 





y = sin A) (100) 
We construct a differential equation as follows 
5 = hos h(t) (101) 
3 = hos h(t) — (4)* sin AC) (102) 
combining equations we obtain 
5 - G/hjp +)?» =0 (103) 


For the special case y = sin wt we obtain ji’ + wy = 0 
and the analog circuit is shown in Fig. 47. 


-ref 0 
cos wt 





w 


Fig. 47—Sine and cosine functions can be generated from 
the solution of a differential equation, 


If h(t) is a function of a variable other then time the 
circuit is somewhat more complex. For example if 
h(t) = x, then we must generate y = sin x, where x and 
its derivatives with respect to time are available. We must 
solve the equation 


a& 
qed = 0 


(104) 
A computer circuit for this function contains multipliers 
in place of potentiometers as shown in Fig. 48. 

When generating integer powers of t, multipliers may 
be used. However, a very convenient method involves in- 
tegrators in cascade with the initial condition on each 
integrator set at zero, These are shown in Fig. 49. 
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Fig. 48—Multipliers are used instead of potentiometers if the 
sine and cosine are to have variable frequency. 
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Fig. 49—Integrators in cascade can be used to generate 
power functions of time. 

Fractional powers of functions can readily be generated 
if the function is not complex or imaginary and if multi- 
pliers and logarithmic elements are available. For example 
the square root, cube root and 2/3 power of x may be 
generated by the circuits shown in Figs. 50, 51 and 52. 


V2 





Fig. 50—This circuit generates square roots if the value of x 
are positive. 
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Fig. 51—This circuit generates the cube root. 





Fig. 52—Fractional exponents also can be generated. 


Variable powers may also readily be generated. The 
solution of y = x* where both x and z are variables can 
be obtained by using logarithmic function generators, The 
function may be generated using the circuit in Fig. 53. 


The logarithm of a function may be generated by the 
use of function generators as well as by the solution of 
differential equation. For example to obtain z = Inx we 
take the time derivative of both sides, 


dz/dt = (1/x)dx/dt 
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Fig. 53—Logarithmic function generators are used for variable 
exponent where x and the powers of x are positive. 





Fig. 54—A time derivative can generate a logarithm. 


and obtain the computer circuit given in Fig. 54. 

One must be careful to scale carefully when using the 
differential equation technique of function generation. 
This is especially important when nonlinear equations are 
solved since nonlinear computer elements may become 
overloaded within the range of interest of one of the 
variables. For example, suppose we wished to generate 
z = t¥ by the differential equation technique. This may 
be written as In z = X (1/3) In t, and by differentiating 
we obtain 


(1/z)(dz/dt) = X 1/3¢ 
X 2/3t 


(106) 
(107) 


or 





This cannot be instrumented directly on the computer 
since at the start of the problem t is zero and a divisor 
cannot be allowed to go to zero. However, one can gen- 
erate 


z= (t+) fora >o (108) 


by the differential equation technique. The equation to 
be solved is 


x (109) 





30 +a) 


and the computer circuit is given in Fig. 55. 








3(t+a) 














Fig. 55—Proper scaling must be used when generating a func- 
tion with the differential equation technique. 
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Function Generation Devices. In some types of prob- 
lems the engineer may have to generate functions of the 
form y = f(x) which may be discontinuous, empirical, 
or nonanalytic and thus not readily obtained as a solu- 
tion of a differential equation. There are several types of 
devices which have been developed to produce the desired 
function or an approximation to it. Some of the more 
often used function generators are: 

® Servo curve followers 

© Servo motor driven tapped potentiometers 

© Servo motor driven potentiometers wound with re- 

sistance wire so that resistance is a particular function 
of shaft position 

® Diode function generators which approximate func- 

tions with straight line segments. 

The tapped servo motor driven potentiometer and the 
diode function generator are the most commonly used 
function generators. Both of these types of function gen- 
erators approximate the function of interest by a series of 
straight line approximations. 


A tapped potentiometer which may be used to approxi- 
mate a function f(x) by straight line segments is shown 
in Fig. 56. Voltages are impressed at the taps and when 
the slides are in some position x between two taps the volt- 
age f(x) is a linear interpolation of the voltages between 
the two taps. Servo motor driven tapped potentiometers 
are usually provided with from 10 to 20 taps so that a 
function may be approximated with from 10 to 20 straight 
lines. Fig. 56 shows how the function represented by the 
solid curve may be approximated by the dotted straight 
lines. 


f(x) 


Servo 
Motor 
Drive 





Fig. 56—Tapped potentiometers can approximate a function 
by a series of linear functions. 


The taps on most servo motor driven potentiometers 
are usually located a fixed distance apart. This feature 
may limit the usefulness of this device for a complex func- 
tion which has several points of inflection which are close 
together. 

Servo motor driven potentiometers wound with resist- 
ance wire so that resistance is a function of the potentiom- 
eter shaft position are usually found only in analog systems 
that require repeated use of a particular complex func- 
tion. This type of device can not be varied and must be 
carefully chosen and scaled for the particular application. 

Modified X-Y recorders are often used as a curve fol- 
lower type of function generator. The function of inter- 
est is plotted using conducting ink and by using a special 
signal pickup on the X-Y recorder the functional relation 
may be retraced and the resulting signal used in the ana- 
log computation. 

The use of function generating devices employing 
servo motors requires careful time scaling of the analog 
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Fig. 57—A diode function generator is faster than servo- 
motor-driven potentiometers. 
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Fig. 58—The simplest diode function generators do not have 
a means of varying the breakpoints. 


circuit so that the problem is scaled within the frequency 
limitations of the servo motors. This means that fast re- 
petitive operation employing oscilloscope display can not 
be used, and the solution must usually be recorded with 
a servo recorder. 


Diode Function Generation. Much of the servo motor 
driven function generation equipment has been replaced 
by diode function generators in modern analog computers. 
Although diode function generators represent functions 
by straight line segments just as servo motor driven po- 
tentiometer function generators, the diode function gen- 
erator is an all electronic system and is therefore much 
faster than servo motor driven systems. 

These are two general types of diode function genera- 
tors: fixed and variable breakpoint generators. The volt- 
age at which two line segment approximations intersect 
(a breakpoint) is fixed by the computer manufacturer in 
a fixed breakpoint function generator, while the engineer 
solving the problem may decide, within certain limita- 
tions, the location of breakpoints for a variable breakpoint 
function generator. The diode has such an important 
function in diode function generation, that it is essential 
its operation be understood. In a previous article in this 
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Fig. 59—Complex lines can be simulated by adding the out 
put from several diode function generators. 
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Fig. 60—The relation between viscosity and shear rate is di- 
vided into segments. 


series it was shown that if the plate voltage of a diode is 
larger than the cathode voltage the diode permits a cur- 
rent to pass through it. If the reverse is true, the diode 
acts as an infinite resistance and no current can pass 
through it. Fig. 57 illustrates the properties of a diode. 

In the circuit of Fig. 57 the diode acts as a switch. 
When the applied voltage ¢, is less than the bias voltage 
€o, the diode acts as an infinite resistance. At values of 
applied voltage equal to or greater than the bias voltage 
the diode acts as if it has zero resistance and current may 
flow. The entire voltage drop 2 due to the flow of current 
though the diode appears across the resistance R as 

@ = (4 — 4)/R (110) 

The solid curve in Fig. 57 illustrates how e, would vary 
with an increase in e, if the bias voltage were set at ey for 
an ideal diode. Actual diodes have characteristics more 
closely represented by the dashed curve. The value of e» 
of bias voltage in Fig. 57 is called a breakpoint. The grad- 
ual rather than sharp change of slope at the breakpoint 
of the diode characteristic curve permits a more accurate 
representation of many functions which do not have sharp 
changes in slope. 
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Function Simulation. The process of simulating a func- 
tion with diode function generators is one which involves 
setting the slopes and breakpoints of the various line seg- 
ments which make up the approximation. It is clear that 
by setting e, in Equation 110 a breakpoint may be fixed, 
and by setting R a slope may be fixed. The simplest diode 
function generators do not have a means of varying the 
location of breakpoints. Fig. 58 shows a simplified fixed 
breakpoint diode function generator. 

A 10-volt computer with a 10-diode fixed breakpoint 
function generator would normally have breakpoints at 
I-volt intervals. The slopes of each of the line segments 
are set by adjusting the values of resistors in the circuits 
shown in Fig. 58. Fig. 59 shows how the slopes of the first 
four line segments might be set to obtain the composite 
function shown as a dashed line. 

It should be noted that as the independent variable x 
increases from x, the dependent variable f(x) follows the 
curve starting at x, until x reaches the value x,. At this 
point the next diode circuit conducts and the value of 
f(x) is the sum of the curve starting at x, and the curve 
at x,. As x continues to increase each line segment adds 
in as the value of the breakpoint voltage corresponding 
to the function is reached. Most diode function generators 
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Fig. 61—This circuit solves the example problem. 


contain an adjustment which permits the operator to 
translate the entire function f(x) between the two ex- 
tremes of machine voltage. 

A fixed breakpoint diode function generator should only 
be used for functions that do not have sharp or closely 
spaced changes in slope. The variable breakpoint diode 
function generator has greater utility than the fixed break- 
point type because of the ease with which breakpoints may 
be adjusted to conform to the variations of the function. 
The procedure for the adjustment of breakpoints and 
slopes may vary depending upon the construction of the 
computer. Before attempting to generate a function by 
this technique the computer manual should be consulted 
for procedural details. 

A simple method of setting up and checking a function 
is to follow these steps: 

Step 1. Plot the desired function on a sheet of graph 
paper to a large scale. 

Step 2. Locate breakpoints. If a fixed breakpoint func- 
tion generator is used, the breakpoints are usually located 
at even voltages. If a variable breakpoint function gener- 
ator is used, convenient points to locate breakpoints are 
at points of inflection. A great deal of judgment is needed 
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Fig. 62—How viscosity changes with startup time. 
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Fig. 63—Power requirements increase during startup. 


TABLE 12—Diode Function Generator Settings 
—_—.] 


rs 


1.000 
0.998 
0.977 
0.932 
0.865 
0.780 
0.667 
0.537 
0.462 
0.412 
0.378 


ey 


2. 
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TABLE 13—Pot Settings 
ee 





Pot Setting 
oO 0.788 
5 0.884 
6 0.176 
35 0.011 
50 0.021 
51 0.834 





to determine the best location of breakpoints. It should 
be noted that breakpoints must be separated by a specified 
voltage depending upon machine construction. Otherwise 
the diodes will interact. 
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Step 3. A table should be prepared listing breakpoint 
voltage and the value of the function at each breakpoint. 
It is some times helpful to approximate the function on 
the graph with straight lines between the breakpoints to 
obtain an idea of the shape of the function to be gen- 
erated. 

Step 4. Potentiometers may now be adjusted to obtain 
the proper value of f(x) corresponding to each break- 
point voltage x. This may be conveniently done in the 
operate mode. 

Step 5. A time ramp may now be applied as the inde- 
pendent variable and the generated function plotted with 
an x-y plotter. The plotted curve may be compared with 
the desired function and potentiometers adjusted until 
the two curves conform as well as possible. 


An Example. Now consider the use of a diode function 
generator in a computer for a problem in which a non- 
Newtonian fluid is to be mixed. We are interested in the 
power required to mix this fluid and in particular the 
power required as the mixer starts and accelerates to its 
final speed. We assume that the mixer speed characteris- 
tics can be approximated by an exponential. The relation- 
ship between viscosity » and shear rate y for the non- 
Newtonian fluid is shown in Fig. 60. This relationship 
will be approximated with a 10-segment fixed breakpoint 
diode function generator. We will assume that there is a 
linear relationship between mixer speed N and shear rate 
y of the form N = ky. We assume that the relationship 


\a /y72 y\d 
Pape ela) (as 2 
pN'pE K a oR 4 (111) 
describes the mixing process, where: 
P = power 
JN = speed 
D = impeller diameter 
p = fluid density 
# = viscosity 
a,b = negative constants 
XK = constant 
This relation may be reduced to 
P=K' Nyt (112) 
where Ri = Kp pet 4/3 (113) 
¢ = a+2b+3 
a=-¢ 


We may solve this equation with the computer circuit 
shown in Fig. 61. 

By using Fig. 60 we may determine the functional 
values at each breakpoint. The function generator used 
in solving this example problem is of the fixed breakpoint 
type with breakpoints at 10 equally spaced intervals. The 
values of the function at each breakpoint are contained 
in Table 12. The pot settings used in this simulation are 
given in Table 13. The results of the simulation shown 
in Figs. 62 and 63 show that the viscosity decreases rapidly 
to a steady value and the power required increases al- 
most exponentially. 


Indexing Terms: Analogs-9, Calculus-6, Circuits-10, Computations-4, Com- 
iters-9, pao Diodes-10, Electricity-10, Engineering-4, Functions-6, 
rograming-7,10, Simulation-4. 
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PART 8: MEMORY AND LOGIC-Only 
modest increases in these will give an 
analog computer a marked increase 
in problem-solving capabilities 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


The memory and logic aspects of analog computer 
implementation are of particular significance to the prac- 
ticing engineer because of the marked increase in prob- 
lem-solving capacity which results with even a modest 
increase in memory and logic capability. Trial and 
error procedures may be automated, the quantitative ef- 
fect of parameter changes may be automatically deter- 
mined, and operating conditions may be changed auto- 
matically during a simulation according to criteria 
specified by the programmer. With such options, among 
others, capable of being exercised during normal or re- 
petitive analog operation, the analog computer becomes 
a very rapid and powerful tool for the solution of many 
problems of particular significance to chemical and pe- 
troleum engineers. 

This particular area of analog computation is quite ex- 
tensive and, with increasing sophistication, merges into the 
large field of techniques better described by the word 
hybrid. However, the material which will be presented 
here shall be restricted to the discussion and illustration 
of the use of analog memory and logic capabilities which 
are now commercially available on moderate size analog 
facilities. 


The Distinction Between Analog Signal and Logical 
Signal. Before discussing the individual components used 
to achieve high speed analog memory and logic, it is 
convenient to first make the distinction between an analog 
signal and a logical signal. Thus far in this series, dis- 
cussions have revolved around analog signals. These are 
the outputs from common analog components such as the 
summer, integrator, multiplier, etc. and have values which 
may range from plus reference to minus reference. Logi- 
cal signals, on the other hand, have only two possible 
voltage values. For example, either 0 or 5 volts. These 
two possible values may be referred to in several equiv- 
alent manners as follows: 


© HIGH and LOW respectively if one prefers to refer 
to the relative voltage levels. 
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© Binary 1 and Binary 0 respectively if one prefers to 
regard the variable value as a binary number. 

© TRUE and FALSE respectively, thereby, retaining 
perhaps more closely the concept of a logical signal. 

In this article the last manner of interpretation will be 
used exclusively. 

The distinction between the two classes of signals is 
necessary because the components on analog computers 
incorporating high speed analog memory and logic fall 
into four main categories: 

© Components which operate with only analog signals. 
These have been considered exclusively to date in this 
series. 

© Components providing the interfacing from analog 
signals to logical signals. 

© Components which operate with only logical signals. 

© Components providing the interfacing from logical 
signals to analog signals. 

An example perhaps best illustrates the types of com- 
ponents which falls into each of these four categories. 


An Example. Consider the physical system indicated 
schematically in Fig. 64 in which the overflow from Tank 
1 forms the input stream into Tank 2. 


F, cult /min 
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TANK 2 
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Fig. 64—The overflow from Tank 1 is the input stream into 
Tank 2. 


Balances around Tank 1 and Tank 2 yield 


Ax (dhy/dt) = F — Pi — Por (114) 


and 
As (dhn/dt) = Por — P2 (115) 
Assuming P, = f:(/i) and Pz = f2(hz) such that the 
mathematical model will not predict negative values of 
h, or hy nor permit the overflowing of Tank 2, the model 
is complete once the proper conditions are placed on 
Por = Overflow. These conditions are 


Por = F — Prify > handif F> Py (116) 


otherwise 


Por = 0 (117) 
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The simulation of this system using components from 
each of the four component categories resolves into the 
following: 

© Designing an analog circuit assuming the proper 
value of Po, is available. 

© Obtaining a logical signal (L:) which is TRUE if 
hy > h and FALSE otherwise, and a logical signal (Iz) 
which is TRUE if F > P, and FALSE otherwise. 

© Obtaining a logical signal (L;) which is TRUE if 
L, = TRUE and L, = TRUE but FALSE otherwise. 

© Choosing between F — P, and 0 for the value of Por 
depending on the value of Ls. 

A preliminary circuit for the first step, assuming F is 
available, is completed as shown in Fig. 65. 


2 vA, 


a he 

















Way 

















A 





= tm) 











-Por. 











P= talhed 








Fig. 65—This preliminary circuit determines the relation be- 
tween feed rate and the heights in the two tanks. 


The logical signals L, and L, may be obtained using 
a component, known as the electronic comparator, hav- 
ing the symbol and characteristics given in Fig. 66. The 
indicated characteristics for this and other components 
presented may vary from machine to machine. It is sug- 
gested that compatability be established before using the 
circuits presented. The signals a, and a, are input analog 
signals, and L and L are output logical variables. 
































a +g >€| TRUE | False 
a+ 0, <€| False | TRUE 


Fig. 66—The logic signal comes from an electronic compar- 
ator. 





In Fig. 66, ¢ is the small voltage difference required for 
activation. It is of the order of a millivolt. The symbol L 
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is called the logical complement of L and is obtained 
from an “Inversion” of a logical signal. 

It is worth noting that the bang-bang circuit given in 
Part 6 of this series can be used as a comparator by 
making the proper choice for the lower and upper volt- 
age levels. 


Using an electronic comparator, the desired signal Ly 
can be obtained by letting a, = hy, a. = —h,andL = Ly. 
Similarly L, can be obtained by letting a: = —F, a; = Py, 
and L = Toes 

The logical signal L, is obtained by operation on L; 
and Z, using a component known as the AND gate. The 
symbol and characteristics are given in Fig. 67 and /, and 


J, are input logical signals and L and L are output logical 


signals. L is referred to as the AND signal and L as the 
NAND (not AND) signal. 
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FALSE FALSE racse | TRUE 








Fig. 67—The AND gate gives the logic signal Ls. 


Letting 1, = Ly, lk = Zp, and L = Ls, the third step of 
the example is completed. 


The final step incorporates the use of a component 
known as the electronic switch. Its characteristics and 
symbol are given in Fig, 68 where L is an input logical 
signal, a; is an input analog signal, and a, is an output 
analog signal. 


" 1 LIS THEN 
9 oo | TRUE a9 = 9) 


FALSE | UTPUT 1S GROUNDED 

















Fig. 68—An electronic switch is used to let the input logical 
signal control the analog output signal. 


As noted in Part 6 of this series, relay switches oper- 
ating with only analog signals are available. A schematic 
diagram, together with operational characteristics, of one 
such switch is given in Fig. 69. 
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Fig. 69—Analog signal switching can also be accomplished 
with a relay circuit. 
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Fig. 70—This electronic circuit is equivalent to the relay sig- 
nal switching shown in Fig. 69. 


The electronic equipment required to achieve the 
equivalent is given in Fig. 70. 

For normal speed of analog operation the relay switch 
is usually quite satisfactory, For repetitive, high speed op- 
eration, however, the wear and required switching time 
can be excessive. The electronic switch has no such limi- 
tations and its use is highly recommended over the use 
of the relay switch for those cases in which a rapid 
switching rate may be required. 


The results of the four steps may now be combined 
in order to obtain a circuit for the example problem. 
‘The circuit, exclusive of scaling, is given in Fig. 71. 
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Fig. 71—This combined circuit for the overflow program does 
not include scaling. 


The reader may verify that the example could also 
be solved using a relay switch and a limiter in order to 
impose the conditions on Por. Speed limitations and diode 
rounding would, however, restrict accuracy of solution, 
particularly for high speed, repetitive operation. 


The OR Operations 

In the same manner in which components operating 
with only analog signals can be connected in a variety 
of ways in order to achieve a desired operation, so may 
the components which operate with only logical signals. 
Two examples are the use of AND gates and logical in- 
verters in order to achieve the OR operation (Fig. 72) 
and the exclusive OR operation (Fig. 73). The validity 
of the circuits may be readily verified by sequentially 
following the possible sets of input values through the 
circuit. 
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Fig. 72—The symbol on the right is a simplified schematic 
of the one on the left. This is an OR operation. 
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Fig. 73—L, is the result of the exclusive OR operation on the 
I, and |; signals. 


Point Storage of Signals. In addition to the electronic 
comparator, the AND gate, and the electronic switch, 
there is a fourth component which is essential for a basic, 
high speed memory and logical capability. This compo- 
nent is the track and store unit. 

The track and store unit can be classified as being a 
component which provides the interfacing from logical 
to analog signals. More specifically, it permits a point 
value of an analog signal to be stored at a specific time 
or condition provided that a logical signal change corre- 
sponding to that specific time or condition can be syn- 
thesized. The symbol and characteristics are detailed in 
Fig. 74 where a; is the input analog signal, a, is the out- 
put analog signal, and L is the input logical signal. Dur- 
ing the condition referred to as track, an internal capac- 
itor is charged by the input analog signal. During the 
condition referred to as store, the input analog signal is 
not connected to the capacitor and the charge which ex- 
isted on the capacitor the last time the input was con- 
nected provides the output analog signal. 
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Fig. 74—This is a track and store unit. 


Track and Store Units. There are numerous ways in 
which the track and store units can be used to great 
advantage during the solution of engineering problems on 
the analog. Several of these are outlined below. 


Use of a Single Track and Store Unit. It is frequently 
the case during the completion of process designs that a 
particular value of a particular analog variable is of pri- 
mary concern. This may be the maximum value of the 
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femperature in a reactor, the length of a reactor required 
to achieve a desired conversion, the percent recovery ob- 
tained with a given gas absorber, etc. The process of 
finding this particular value can be automated using a 
single track and store unit provided the proper logical 
signal can be synthesized. The restriction is that the logi- 
cal signal must change only at the desired condition. 

Suppose, for example, that it is desired to determine 
the maximum temperature in a flow reactor. First an 
analog circuit is designed for calculating T versus reactor 
bed length, x. If preliminary calculations have shown that 
the temperature increases with length to the maximum 
and then continuously decreases, the maximum tempera- 
ture may be automatically determined as shown in Fig. 75. 

During a single computer run, A will be equal to T 
until dT/dx = 0 and thereafter A will remain constant 
BUT was. 
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Fig. 75—A circuit to determine the maximum temperature of 
@ flow reactor bed. 


Transfer of a Point Value from One Analog Circuit 
to Another. The output of the track and store unit can 
be used for further analog calculation in the same man- 
ner that the output from a summer, integrator, etc. may. 
Consider, for example, the case in which a reactor- 
absorber system such as shown in Fig. 76 is to be sim- 
ulated. 
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Fig. 76—The output composition y, depends upon all of the 
other compositions and dimensions shown here. 


Given an analog circuit for finding y vs. x, by integrat- 
ing dy/dx, for the reactor and another for finding ¥ vs. z, 
by integrating dy/dz, within the absorber, assume that it 
is desired to determine the absorber height required to 
reduce yz to f(y:), 0 <f <1, when y, and X are known. 
The process of first solving the reactor and then the ab- 
sorber is to be automated by making one run to deter- 
mine y|--x from the reactor (during which time the ab- 
sorber results are ignored) and then a second run in 
which the absorber equations are to be solved for zZ. 
Determining y|,-x is indicated in Fig. 77. However, upon 
placing the computer in reset mode so that the absorber 
circuit may be initialized before solving, x drops to zero 
and the track and store unit begins to track y = yo. Asa 
consequence, y|,-x is lost. In order to achieve the transfer, 
two track and store units must be used in series as shown 
in Fig. 78. This arrangement is commonly referred to as 
the “bucket bridgade.” 

The sequence of events which occur for the example 


REPRINTED FROM HYDROCARBON PROCESSING 






































REACTOR ne atlz=xj tat x, 
Yo circuit 
-«—_ [+ 
x Aen SS 


ore 


Fig. 77—The reactor-absorber in Fig. 76 is simulated. 
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problem during the process of making two runs are as 
follows: 

® During the first run, y’ is initially equal to —y since 
Lis TRUE. y” is equal to the value of —y’ the last time 
L was TRUE since L is now FALSE (This value is zero 
in this case if the computer was placed in the pot set, 
the reset, and then the operate mode when the run was 
started). 

© When x > X occurs, Z becomes FALSE and f= 
—y|z=x. However, L becomes TRUE and y” represents 
the tracking of —y’. Consequently, y” = y|.a-x. 

© Upon placing the computer in the reset mode in 
order to begin the second run, ZL becomes TRUE and 
x = —yo. L is now, however, FALSE and thus y” is equal 
to the value of y’ the last time L was TRUE. (y”” = yleax 
obtained during the first run.) 

© When the computer is placed in the operate mode 
for the second run, y” will remain = (y|,-x obtained dur- 
ing the first run) until x > X, at which time y” = (y| eax 
obtained during the second run). 

Using the bucket bridgade to insure the proper trans- 
fer of y|z-x, the reactor-absorber system may be solved 
automatically in two runs by incorporating the circuit in- 
dicated in Fig, 79. The signal A will be equal to the de- 
sired value of Z during the latter portion of the second 
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Fig. aa complete simulation of the reactor-absorber shown 
in Fig. 76. 


Iterative Solutions Using Point Analog Values. In the 
above example, the bucket bridgade was used to transfer 
a point analog value from one analog circuit to another 
during the reset period between runs. A closely analogous 
situation arises in those cases in which point analog values 
are used in iterative solutions, 

Consider as an example the calculation of the percent 
recovery in a gas absorber of known height. (Refer to 
Fig. 76 for a schematic of an absorber). Assuming that 
a differential equation of the form dy/dz = f(x,y) is 
available for the absorber, the value of y, and x, must 
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both be known in order to integrate the differential equa- 
tions by letting z be the computer operating time. In the 
general case y; and x2 are known but x, must be calcu- 
lated by an overall balance assuming a value for ys. 
That is: 


x1 = x2 + (V/L) (1 — 22) (118) 


for the simplified case in which V/L is a constant. The 
assumed value of 92 must then be checked and the cal- 
culation repeated using another assumed value if a sig- 
nificant error is observed. In many cases, a process of 
iterative substitution can be used for this type of calcu- 
lation. The manner in which a bucket bridgade can be 
used for the iterative solution of the recovery in an 
absorber is indicated in Fig. 80. 
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Fig. 80—A iterative solution can be employed. 





The signal A, during the period of calculation in which 
z < Z, will be equal to the value of yz which was calcu- 
lated during the last run. When z > Z, ye is updated. By 
then placing the computer in reset, a new run may be 
begun using the updated value of y2. Such an iterative 
calculation can be completed automatically by placing 
the computer in repetitive operation mode during which 
time the computer is automatically cycled between the 
reset and the operate modes. 


Single Parameter Sweep. Very seldom is the practicing 
engineer’s primary concern the solution of a neatly de- 
fined mathematical model. In the general case, assump- 
tions of dubious physical validity must be investigated 
and engineering judgment exercised in the selection of 
parameter values. Consequently, an engineer may be 
called upon to propose several models and, for each 
model, to investigate quantitatively the effect of param- 
eter changes. In this section, a technique is discussed 
which permits the automation of the steps required to 
carry out a single-parameter sweep. 

Suppose that an analog circuit for the system under 
consideration is available and that a bucket bridgade has 
been designed so that a desired point value of an im- 
portant analog signal can be retained for plotting versus 
the parameter P. (See Fig. 81) 


L 
p——J anatoc 1 
CIRCUIT 


Assuming that the analog circuit is in the repetitive op- 
eration mode with an operate period of at least Z, which 
is in turn greater than the time required for L to change 
to FALSE, it is desired to automate the following pro- 
cedure: 


Fig. 81—A bucket bri- 
gade for a single param- 
eter. 
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© Maintain P at a fixed value during the operating 
period for Z time units. 

© When the computer resets, increment (or decre- 
ment) the fixed value of P by an amount AP. The value 
of A will then be plotted versus P while the computer is 
in repetitive operation. 

The solution utilizes a configuration of track and store 
units called the accumulator which is given in Fig. 82. 
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Fig. 82—The accumula- 
tor circuit. 




















TABLE 14—The Sequence of Events in an 
Accumulator 






Event 





Assuming L is initially TRUE and a, = Py, the se- 
quence of events occurs as L changes as shown in Table 
14. Thus from the time at which L becomes FALSE to 
the time it becomes FALSE again, a, remains constant, 
but az is incremented by AP each time L becomes FALSE. 

In order to automate the steps for the parameter sweep, 
one needs only use the output of the accumulator for 
the parameter value and use the operating time z to syn- 
thesize the logical signal L. The completed circuit is given 
in Fig. 83. 

"ANALOG 


Fig. 83—The single parameter sweep can be automated. 















































Fig. 84—The results from the circuit , 
of Fig. 83. 





P 


A plot of A versus P would have the characteristics 
given in Fig. 84, The desired point values of the analog 
variable as a function of P are given by the points b, d, 
etc. If P is sufficiently small, the plot will appear to be 
continuous. 


Two Parameter Sweep. Consider next the case indi- 
cated schematically in Fig. 85 for which it is desired to 
automatically determine A as a function of P;, in incre- 
ments of AP,, and Pz, in increments AP, by operating the 
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Fig. 85—A two parameter sweep circuit. 
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Fig. 86—The results from the 
circuit of Fig. 85. 


circuit in the repetitive operation mode. It is desired to 
design a circuit which will change P, and P, in the man- 
ner shown in Fig. 86. 

With AP, sufficiently small and AP, rather large, a plot 
of A versus P; will appear continuous with P, being a 
parameter on the plot. 

The circuit in Fig. 87 permits the desired changes to 
be made in P,. 
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Fig. 87—A circuit to determine the change in P,. 


The logical signal L, determines the rate at which P, 
is incremented. The logical signal L, establishes whether 
AP, is positive or negative. If AP, is positive, P; is com- 
pared with P,”; otherwise with P,’. The logical signal Ls 
permits the sign of AP; to be changed. If AP, is initially 
positive and P, — P;” < e, a positive AP, is retained. 
However, when P, — P,” > e occurs, a negative AP; 
is specified. The negative 4 P; causes P; to be compared 
with P,’, The sign of AP, becomes positive again when 
P, — P,’ <e first occurs and P; is again compared with 
P,”. In this manner, the cycling of P, is insured. 

In order to increment P, in the desired manner, the 
circuit given in Fig. 88 has been used. The signals P, 
and Ly are obtained from the circuit in Fig. 87. 


(Pi + PAN Pi 


Fig. 88—A circuit to de- 
termine the change in P;. 
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It is to be noted L, is the result of an exclusive OR 
operation on L, and Ly. Thus L; is TRUE if only L, is 
TRUE or if only L, is TRUE and L; is FALSE other- 
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TABLE 15—The Sequence of Events as P, 
and P, Undergoe Cycling 








Event AP:| Ls | Le | Ls Pa 
+ | False | False | False | Po 
Oe 

(Pi? + Pr)/2 + | False | True | True | Po 
Pi surpasses Pi’ + | True | True | False | P2 + A Po 
A Pi changes sign — | True | True | False] P2 + A Pa 
P1 becomes less than 

(Pr + Pir)/a — | True | False | True | Pa + A Po 
Pi becomeslessthan Pi’ | ~ | False | False | False] Pa + A Pa + A Po 
4 Pr changes sign False | False | False | P2 + A P2 + A Pa. 





wise. The sequence of events given in Table 15 occurs as 
P, undergoes cycling. Assume AP, is initially positive and 
P, is at the minimum. 

The next event will again be the first entry into Table 
15 and thus the desired changes in both P, and P, have 
been logically achieved. 

The final circuit for a two-parameter sweep is obtained 
by combining the circuits shown in Figs. 85, 87 and 88. 


Solution of Boundary Value Problems. Boundary value 
problems can be regarded as the generalization of the 
class of problems mentioned earlier in the section, “Itera- 
tive Solutions Using Point Analog Values.” There an 
unknown input to the analog circuit was identical to a 
desired analog output. In the general case, it is known 
that an unknown input is a function of a desired output 
but the functionality is not available. The manner in 
which the required trial and error calculations may be 
automatically completed is discussed in this section by 








ig. 89—A double-pipe reactor with heat exchange. 


TABLE 16—Nomenclature for the Boundary 
Value Problem 








Variable Definition Units 

A surface area per unit length 8. ft./ft. 
Au crosssectional area available 

for flow of hot fluid sa. ft. 
As crosseectional area available 

for flow of cold sq. ft. 
c mass fraction of reactive species Tbe lb 
Sp heat capacity BTU/(b) (PR) 
BD diameter of inner pipe ft 
E activation energy divided by the gas constant | °R, 
Fa mass flow rate of hot flui W/sec 
Fe fuss Gow nie of cold fel Ib/se 
AH | heat of reaction BTU/Ibe 
k reaction rate constant (eu ft, see)-t 
x length of exchanger 
Tu temperature of hot fluid He 
Te temperature of cold fluid °R 
va velocity of hot flui ft/sec 
ve Velocity of cold fluid ft/sec 
U overall heat transfer coefficient BTU/(sec) 

pe) CR) 

x length along exchanger 
y reaction rate (eu, ft see)! 
° fluid density Tb/eu. ft. 





Special Symbols Used With Variables 

* superscript to indicate expected maximum value of a variable. 

~ overscore to indicate an initial or a desired value. 
engl@t Over variable to indicate the first derivative with respect to exchanger 
length. 
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TABLE 17—System‘s Equations 


nS 
Cold Fluid: 

Te = UA (Te — Tu)/\Fe co) 

where U, A, Fe: cp and Te at + = X are assumed to be known. 


Reaction Rate: 


-zIT 
yoke oH 


or 
y= ETH y/TH) 
where E and k are assumed to be known. 





Reactive Specie 
é=co/viH 


where yy and cat x = are assumed to be known. 





Hot Flutd: 
Py = UA (Te — TH) (Fu cp) + SH 9 c/(vH op) 
where Fy, AH, and Ty at x = O are assumed to be known. 
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Fig. 90—The preliminary circuit for the reactor of Fig. 89. 


considering the simulation of an exothermic reaction in 
a double-pipe heat exchanger. A schematic of the chosen. 
system is given in Fig. 89. Table 16 summarizes the 
notation used. 

If it is assumed that component c undergoes an exo- 
thermic, irreversible, first-order, temperature dependent re- 
action as it flows through the reactor, the differential 
equations presented in Table 17 may be derived in order 
to describe the system. The preliminary analog diagram 
given in Fig. 90 may be designed. 

Since the value of T, and ¢ at x = 0 are assumed 
known, Ty = Ty at x = 0, ¢, and y are known. However, 
the value of T; is known at x = X. In order to solve the 
problem it is necessary to assume a value of T, at x = 0, 
check to see if Te = T- at x = X, and to repeat the cal- 
culation using a new value if an appreciable error exists. 
The relationship between T,, at x = 0 and at x = X is not 
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explicitly known, but from physical considerations it is 
noted that if T, at x = 0 increases, T, at x X will also. 

Thus if T, at x = X is larger than T., Te at x =0 
should be decreased. Similarly if T, at x = X is less than 
T., T; at x = 0 should be increased. It is also noted that 
as the error decreases, the required correction decreases. 
As a first approximation assume the correction to be 
made is proportional to the error. With this assumption, 
letting K be the proportionality constant, it may be veri- 
fied that the addition of an accumulator to the prelimi- 
nary circuit will permit the proper corrections to be made 
as the analog circuit is cycled between the reset and oper- 
ate modes. The resulting unscaled circuit is given in 
Fig. 91. 
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Fig. 91—A basic circuit for solving for boundary values. 
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Fig. 92—A scaled analog circuit gives a numerical illustration 
of the solution of the reactor of Fig. 89. 


Fig. 92 presents a scaled analog diagram which was 
prepared for numerical illustration of this technique. 
Pots Gy, Gs, Gs, and G, have been used in order to scale 
the outputs from the divider, squarer, and multipliers so 
that errors are minimized in the serial use of these non- 
linear components. The values of the required pots are 
summarized in Table 18. 

Using the numerical values given in Table 19, the plot 
given in Fig. 93 was obtained by guessing T, at x = 0 to 
be zero initially, and then cycling the computer between 
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reset and operate. From the plot it is seen that T, — T- 
at x = X was reduced to an insignificant amount after 
only 6 trials. In Fig. 94, the profiles obtained for Ty, Tes 
c, and y after convergence was achieved are summarized. 

The illustrated technique finds general applicability in 
the solution of boundary value problems. It is necessary 
to know only the following: 


TABLE 18—Pot Settings for Fig. 92 





Pot Value 
00 -Te#/(Ts — Tu)* 

OL UA (Te — Tu)*/(Fe cp Te* B) 

02 Ta/Ta* 

Ta*/(Ts — Tu)* 

UA (Ts — Ta)*/(Fu cp TH?) 
Tw*/(Tn* B) 

06 ¥/y* 

G1 Such that (G1 Ta*/Tu)mas ~ REF 
G2 Such that (G1 Tu*/Tw)? (y/y*) (1/ 
G3_ Such that (G1 Tn*/Tx)? (v/y*) (1 
52 G3 Ta*/(G1? Ta*? B) 

08 c/c* 

G4 Such that [y c/(y* c* G4)]mas REF 

09 [AH] y# ce G4/(vm cp Ta®) 

10 y* G4/(vm B) 

———— 





mas ce REF 
(Tn/Ta*) (1/G3)}0 





REF 








TABLE 19—Numerical Values of Parameters 
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Fig. 93—The computer was cycled between reset and operate. 
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Fig. 94—Other numerical values from the computer program. 


© Direction in which a correction must be made when 
an error exists. 


© Approximate magnitude of the required correction. 
The latter is generally determined by trial and error on 
the value of K. 


By placing the computer in the repetitive operation 
mode, the required trial and error is generally completed 
so rapidly that the time for a solution is not noticed. Due 
to the high speed, parameter changes may be made and, 
for all practical purposes, the instdntaneous solution is 
continuously shown on the output oscilloscope. 

In the next article of this series, more sophisticated 
applications and capabilities will be qualitatively discussed 
and a summary will be given of this series in retrospect 
to the developments in the field of digital simulation. 


Indexing Terms: Analogs-9, Circuits-10, Computations-4, Computers-9, De- 
scriptions-8, Electricity-10, " Engineering-4, Logic-9,Memories-9, Program- 
ing-10, Simulation-4. 
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PART 9: FUTURE ROLE— 

The capabilities of the high-speed 
analog computer are compared with 
those of the digital and the 

hybrid computers 


Theodore W. Cadman and Theodore G. Smith 
University of Maryland, College Park, Md. 


THROUGHOUT THIS SERIES, attention was drawn re- 
peatedly to the use of the analog computer for solving sets 
of ordinary—linear or nonlinear—differential equations 
such as those frequently used to describe the dynamic or 
steady-state characteristics of processing systems. The ap- 
plicability of the analog was noted to result from the 
availability of the hardware integrator which permits 
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integration over a wide range of time scales, the con- 
tinuous and parallel nature of analog operation, and the 
availability of analog components which perform other re- 
quired mathematical operations on the voltages used to 
represent the independent variables. These characteristics, 
coupled with provisions for oscilliscope and x-y plotter 
displays made the analog computer an excellent tool for 
the rapid solution and convenient examination of the 
solutions of differential equations. 

The achievement of the full potentials of the analog 
did not, however, immediately follow. Due to the con- 
tinuous and parallel nature of operation, analog com- 
ponents had to be available in sufficient numbers and 
variety in order to simultaneously perform the mathe- 
matical operations required in the solution procedure. An 
analog computer with 10 integrators could not be used 
to solve a problem requiring 11 simultaneous integrations. 
Moreover, an analog circuit had to be designed, patched, 
amplitude and time scaled, and checked out in order to 
insure satisfactory operation. Thus, in simplest terms, hard- 
ware limitations and problem preparation procedures 
were detrimental attributes of the analog. 


Digital Computer Uses. Due to the disadvantages of the 
analog computer, the engineer is naturally interested in 
the evaluation of the possible alternatives before deciding 
to use the analog. In this evaluation, the recent develop- 
ments in digital implementation can be of particular sig- 
nificance. 

Since about the mid 1950s, the engineer has had at his 
disposal programing languages such as FORTRAN, 
MAD, ALGOL, etc., which, together with previously de- 
veloped numerical techniques for integration, differen- 
tiation, interpolation, and so forth, permitted the digital 
computer to be used for the discrete, sequential type of 
numerical calculations with which the engineer was al- 
ready familiar. Such languages permit digital programing 


ee 


An engineer should evaluate alternatives when deciding 


fo use the analog computer 


to be achieved fairly easily and allow the programer 
to make full use of the flexibility of the digital computer. 


Programing the digital computer is, however, quite de- 
tailed and care must be taken to insure the proper logical 
sequence of operations required for solution. Moreover, 
a program, once developed, tends to be specific in nature 
and not readily subject to extensive modification. If one 
neglects the comparison of the time required for solution 
and chooses a problem which is amenable to both analog 
and digital solution, the balance between the time spent 
on analog setup and the time devoted to the development 
of a digital program using a procedure oriented language 
heavily depends on the past experience of the poten- 
tial user. 

The state of affairs with regard to ease of digital pro- 
graming has, however, been markedly changing during 
the past few years. First, there has been a proliferation of 
user-oriented languages aimed at reducing to a minimum 
the complexity of programing certain specific classes of 
problems. Particularly to be noted in conjunction with 
this series are a group of languages (MIDAS, MIMIC, 
DSL/90, just to name three of several dozen) referred 
to as digital simulation languages. Some of their im- 
portant characteristics are: 


© They are specifically oriented towards the numerical 
solution of ordinary differential equations and are prin- 
cipally aimed at simulating analog operation on the digi- 
tal exclusive of the scaling and storage restrictions of the 
analog. 


© Very little previous exposure is required for their 
comprehensive use, 


© The burden of prescribing solution procedure as well 
as the error check of numerical procedures such as inte- 
gration are removed from the user. 


© Their use permits a very significant reduction in the 
time spent in problem preparation over that required for 
a comparable solution on the analog or with a procedure- 
oriented digital language." 





Digital simulation languages are becoming ever in- 
creasingly more popular. Developments will undoubtedly 
continue and potentially can very markedly reduce the 
time required for programing a digital to provide numeri- 
cal solutions to partial as well as ordinary differential 
equations, boundary value problems, and classes of itera- 
tive calculations. The use of one of these languages re- 
quires a large computer because of the storage needed 
for the program. 


The second marked change in digital inplementation 
which can be observed is the continual increasing in the 
degree of user-machine communication. Remote teletype 
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terminals are now commonly used at installations 
equipped with a time sharing digital facility. These sys- 
tems are very convenient for digital “debugging” of 
procedure-oriented based programs and on-the-spot so- 
lution of small to moderate sized problems. The future 
may very well witness the extensive development of user- 
oriented languages for use with teletype terminals, 
thereby significantly enhancing even more user-machine 
communication. 

Other possibilities for the future revolve around active 
work being done to communicate with the digital com- 
puter by means of the written word. Here the objective 
is to develop means by which the user can write his 
problem on specially prepared paper and have the digital 
computer formulate his problem in detailed terms, con- 
struct a program, run the program, and finally deliver 
the results to the user. Communication by means of voice 
is also being seriously considered. 

Of particular interest to design engineers is work 
currently being carried out with the objective of com- 
municating with the digital by means of large oscilli- 
scopes. Technology has probably advanced to the point 
where today it would be possible to develop a system 
which would require the following sequence of steps for 
an engineer to calculate the dynamic response of a plant 
complex: 








© Specify the equipment and interconnections in the 
plant complex, using a light pen, on an oscilliscope. 


© Specify, either on the scope or by means of cards 
and/or tapes, the parameters required, the conditions to 
be imposed, and the results desired. 


© Press a button and have the digital solve the problem 
and display the desired dynamic response(s) on the 
scope. 

Although such a facility would be inherently complex 
in its makeup, the user could be completely unaware of 
the complexity and be directly concerned only with 
drawing lines using a light pen, specifying operating con- 
ditions, and interpreting results. 


Potential of Analog Computers. With the digital 
simulation languages now available and with the possi- 
bility of reducing communication with the digital com- 
puter to typewritten material, oscilliscope displays, the 
hand written word, or even voice, is the use of the analog 
in engineering to become passé in the forseeable future? 
The answer is an emphatic NO! 


In contrast to an analog computer, a digital computer 
is slow. The digital operates on discrete values in a se- 
quential manner in order to obtain a solution. For a given 
problem and procedure the number of discrete values 





51 








which must be used in order to obtain a satisfactory ap- 
proximation has a fixed lower value. Consequently for a 
given machine, the time required for a satisfactory so- 
lution also has a fixed lower value since the values must 
be operated on sequentially according to the procedure 
being used. For the simulation of inherently continuous 
systems, this required time may be excessive. The pres- 
ence of nonlinearities and higher frequency fluctuation of 
problem values tend to markedly increase the required 
time. The size of the problem to be solved has a direct 
bearing on the required time. 

Analog operation is, however, continuous and parallel 
in nature. Its speed of operation can be chosen by the 
user. A solution may be obtained over a period of several 
minutes or in a few thousandths of a second for a great 
variety of problems. The time required is virthally inde- 
pendent of nonlinearities and frequency fluctuations can 
cover a very wide range. Solution time does not depend 
on the number of equations being simultaneously solved. 
As a consequence of these analog characteristics, the time 
required for a comparable solution on the digital can be 
several orders of magnitude larger than the time required 
for an analog solution. Thus the long term potential of 
the analog lies in its ability to simulate continuous sys- 
tems in a very rapid fashion. This speed potential 
coupled with already developed mechanisms for a high 
degree of user-machine interaction establishes the analog 
as a very powerful engineering tool for the foreseeable 
future despite the present and projected advances in 
digital software for programing ease and user-machine 
communication. 


The more recent developments in the field of analog 
computation have been directed towards the more effi- 
cient utilization of the inherent power of the modern 
electronic analog. In Part 8 of this series, four high- 
speed components which permitted automation in several 
aspects of analog implementation were presented. Par- 
ticularly noted was the role of directly or indirectly 
synthesized logical signals in their use. The electronic 
comparator provided means for generating logical sig- 
nals by comparing analog signals; the AND gate per- 
mitted comparison of logical signals; the electronic 
switch permitted an appropriate analog signal to be 
logically chosen, and the track and store unit permitted 
logical storage of the point value of an analog signal. 
Used collectively they were shown to provide means for 
automating trial and error calculations, the imposition 
of conditions during a computer run, and the process of 
obtaining parameter sweeps. These applications, together 
with the numerous others possible, relieve the user of the 
necessity of making many tedious changes in the analog 
parameters or circuit during or between analog runs and 
permit the effective utilization of analog solutions which 
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may be obtained at a rate of up to about 1,000 per 
second. 

Even the impressive applications which have been dis- 
cussed are, however, overshadowed by the capabilities of 
the larger modern analog computers commercially availa- 
ble. It is our objective here to indicate briefly some of 
these more sophisticated capabilities. 


Automatic Control of Integrator Operation. Essential 
to the more effective utilization of the speed potential of 
the analog is the ability to automatically control its opera- 
tion. We have already seen the additional degree of con- 
trol which the electronic switch permits and the usefulness 
of the track and store unit for controlled memory of point 
values of analog variables. The basic control of an analog 
circuit lies, however, in the integrators. 


On smaller analogs the operation of the integrators 
may be specified manually to be reset, hold, or operate or 
the integrators may be operated in repetitive mode, cycling 
between reset and operate for specified periods of time. 
On the modern larger analogs, integrators may be logically 
controlled in a manner such as that indicated in Fig. 95. 
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Fig. 95—Logical control for an integrator. 


With the logical signals 1G and OP, the user has the 
ability to automatically control the state of the integrator. 
The logical signals 4, and l, permit the speed of integra- 
tion to be varied from a time constant of 1 second, cor- 
responding to normal speed of operation, to 10,000 times 
faster than normal. 

One can use the resulting capabilities in completing a 
parameter sweep, for example. The parameter to be varied 
could be the output from an integrator which remains in 
operate with a time constant of 1 second while the system 
being simulated is run in repetitive mode with a con- 
siderably smaller time constant. 


During one analog run simulating the system, the 
change in the parameter will be negligibly small. During 
several thousand runs (which may take only a few sec- 
onds), the parameter is swept, however, and a parameter 
sweep is completed with less equipment requirements than 
those noted in Part 8 of this series. 

Reversing the method of operating the system and the 
parameter generation permits complex parameter genera- 
tion, such as representing an equilibrium state in the 
system by the steady-state solution of a set of differential 
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equations, to be achieved automatically and continuously. 
The user might also automatically change the time con- 
stant so as to operate very rapidly during a trial and error 
solution and then slow down the solution rate after con- 
vergence so that a permanent record may be made on the 
logically controlled x-y recorders available. The possibility 
of multispeed calculations and automatic selection of the 
mode of individual integrator’s operation thus open a 
whole new realm of analog applications which can be 
used to great advantage in engineering. 


Patchable Digital Logic. A necessary prerequisite for 
the automatic control of integrators, electronic switches, 
and track and store units is the synthesis of the appropri- 
ate logical signals. The electronic comparator provides 
one means for generating logical signals but is dependent 
on an instantaneous comprrison of analog signals. The 
AND gate provides mea:is for continuously comparing 
logical signals. These, together with others to be briefly 
discussed later, constitute what is referred to as the patch- 
able digital logic capability of the modern analog and 
provide the user with the capability of logical signal syn- 
thesis. 


The flip-flop is an essential element in patchable digital 
logic. In order to insure proper logical operations for cir- 
cuits using flip-flops, their operation on many analog 
computers is governed by carefully timed pulses (10° per 
second, 10° per second, etc.). Essentially each pulse can 
be thought of as a spike. During the rising part of the 
spike the flip-flop “decides” what it is to do. During the 
falling part of spike, the action is completed. The sche- 
matic and characteristics of operation of a flip-flop are 
given in Fig. 96. 


SET 

ENABLE- 

RESET- 
———— 


"During Rising Then Following 
Part of Spi Falling Part of 


OUTPUT 















False No change 
True Faise True 

True True False 

True False No change 
True True Changes state 


——— 
Fig. 96—The flip-flop is an essential element. 


As one example of the flip-flop’s applicability suppose 
that an electronic switch is to be thrown the first time that 
an analog variable exceeds a certain value and is to re- 
main thrown throughout the remainder of an analog run. 
(Such an application would arise if there were a change 
in the method of operating a physical system once a cer- 
tain condition was reached.) In this case, the enable 
input on the flip-flop could be made TRUE. Then during 
the IC period of the integrators, the set input would be 
made FALSE and the reset input made TRUE so that 
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the output would be initialized to FALSE. During the 
operate period, the reset input would be made FALSE 
and the set input would be the output from an electronic 
comparator governed by the desired analog comparison. 
The first time the output from the comparator becomes 
TRUE, the flip-flop will set on the falling part of the 
next spike and the output will remain TRUE for the 
remainder of the run regardless of the subsequent result 
of the analog comparison. 


Further insight into the operation of a flip-flop can be 
obtained by examining a schematic representation of a 
component known as the differentiator such as is given 
in Fig. 97. In this case, assume the enable input to the 
flip-flop remains TRUE. 









































SET as 
OUTPUT | 
INPUT | OUTPUT + OUTPUT 3 
OUTPUT 2 
RESET 


Fig. 97—Schematic representation of a differentiator. 


To understand how the differentiator works, assume 
that input and output are initially FALSE. ‘The sequence 
of events which occurs when the input becomes TRUE 
for a length of time and then becomes FALSE again is 
shown in Table 20. 


TABLE 20—Sequence for a Differentiator 






Timing Input | Set 















Rising spike False | False False 
Falling spike True | True True 
Between pulses....} True | True True 
Rising spike True | True True 
Falling spike......| True | True False 
While Input 

remains TRUE... True | True False 
Falling spike False | False True 
Between pulses. False | False True 
Rising spike. . . False | Flase True 
Falling spike False | False False 
While Input 

remains FALSE. .| False | False False 





Examination of the sequence of events shows that 
Output 1 becomes TRUE for only one pulse after the 
input becomes TRUE. Output 1 is termed the leading 
edge differentiation of the input. Output 2 becomes 
TRUE for only one pulse after the input becomes FALSE 
and is termed the trailing edge differentiation of the 
input. Output 3 is the result of the OR operation on 
Output 1 and Output 2 and thus is TRUE for only one 
pulse after the input changes from TRUE to FALSE or 
from FALSE to TRUE. 


The differentiator’s outputs are used when it is desired 
to note the change of logical signals by a TRUE signal 


53 





only one pulse long. They are frequently used in conjunc- 
tion with a component known as the monostable timer. 
The monostable timer is given the symbol shown in Fig. 98. 





AT 


ea a7 ee OUTPUT 


Fig. 98—The symbol for a monostable timer. 











If the input to the monostable timer becomes TRUE, 
even for one pulse, the output from the monostable timer 
will remain TRUE for a period of time specified by A T 
and then become FALSE. The period A T is generally 
manually adjustable from a few microseconds to upwards 
of 100 seconds. The monostable timer thus permits the 
user some timing capability during analog operation. 
One might, for example, use the monostable time to stop 
an analog program while the pen on a plotter is raised 
or lowered, to run a particular portion of an analog circuit 
for an approximate period of time, or to allow time to 
pass so that the initial conditions may be established on 
certain integrators. 


For more precise timing operations than are possible 
with the monostable timer, analogs incorporating a patch- 
able digital logic capability are equipped with counters 
and registers, Counters permit, for example, the genera- 
tion of a TRUE signal only after a logical signal has been 
TRUE for a specified number of pulses. Registers permit 
a sequence of events to be carried out according to a 
rigidly determined time table. 

One of the documented applications of the use of patch- 
able digital logic which has particular applications in 
control studies involves the completion of a double param- 
eter sweep with the automatic plotting of the values of 
the parameters for which the system is stable."* In this 
particular application, the system is tested for stability 
using many thousand possible parameter values in the 
matter of a few minutes. Such an application of modern 
analog techniques can but impress on the engineer the 
power of the modern analog and the very significant time 
savings which are possible with the analog in the simula- 
tion of continuous systems requiring numerous iterative 
and complex calculations. 

Other groups of problems for which patchable digital 
logic is particularly useful are those requiring dynamic 
optimizations, the fitting of multiple boundary values, 
complex sequencing of analog simulations, and in short, 
those problems in which the full potential of the analog 
can not be achieved because of the insufficiency of ordi- 
nary methods for controlling analog operation. These 
classes of problems are among the most complex in engi- 
neering and are more often than not, those for which 
possible alternative methods of solution are very time 
consuming. 


The Hybrid Computer. An analog computer equipped 
with patchable digital logic capabilities is frequently 
called a hybrid computer. In this series, however, the term 
hybrid computer has been reserved to describe the even 
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more powerful computational facility which results when 
an analog computer and a digital computer per se are 
operated as one system, being coupled by a linkage sub- 
system which permits the transfer of operating and control 
information to and from the one computer system to 
the other. 

The need and usefulness of the hybrid are many. From 
the digital standpoint, the inclusion of the analog permits 
very significant reductions to be made in the time which 
would be required to simulate continuous systems and 
systems described by differential equations in particular. 
From the analog standpoint, the digital provides the 
severely needed storage capacity which the analog lacks 
and provides a means for completing arithmetic opera- 
tions which may require an excessive amount of analog 
capacity or a degree of accuracy for which the analog is 
insufficient, Examining the hybrid system as a whole, it 
is observed that the combined use of the analog and 
digital, with appropriate interfacing, permits a system to 
be obtained which retains the advantages of each com- 
puter subsystem while minimizing the deterimental effect 
of their individual disadvantages. 

The achievement of the full advantage of each com- 
puter subsystem places a great deal of importance on the 
ability to transfer, to receive, and to use received infor- 
mation. The linkage is understandably very important and 
its design has been proven to be more complex than may 
be first realized.**** Essential in the analog subsystem 
are means for effective control. Here the capabilities pro- 
vided by patchable digital logic are particularly conve- 
nient. In addition, the availability of servo-set potentiom- 
eters on the analog are required in order to permit the 
full utilization of information available from the digital. 

Particularly important to effective digital utilization is 
the availability of a full complement of software specifi- 
cally designed for use in the hybrid system. These desir 
ble characteristics, as well as others, too numerous to 
mention, are now available on commercially produced 
hybrids.*® 





A few examples are now examined to show the applica- 
bility and application of hybrid computers in chemical 
and petroleum engineering: 

One of the important general uses of the digital in a 
hybrid system is to facilitate problem setup and checking 
out of the analog computer before operation. Satisfactory 
digital software can be used during this initial period to 
complete the calculation of scaled pot settings, set the pots 
by means of servo-set pots on the analog, automatically 
complete a static check and inform the user of the most 
likely trouble spots should one of the provided checks fail. 
Such a capability can reduce programing and setup 
time on the analog by two-thirds.’ 

During the actual solution period the digital may be 
used to contro] the analog, to readjust pot settings after 
an analog run, and to solve a portion of the problem 
itself. Such a capability permits the simulation of complex 
pieces of processing equipment. Fixed bed catalytic re- 
actors“ and tubular reactors'® (which may require simu- 
lation by several hundred continuous stirred tank re- 
actors), multicomponent separation, and a reacting 
distillation column™ have been simulated by using the 
digital for intermediate calculations and control of analog 





LEARN ABOUT ANALOG COMPUTERS . . . 


“subroutines.” The digital’s operation sequentially chooses 
a continuous stirred tank reactor or distillation tray 
for analog simulation, scales the analog circuit, specifies 
the initial conditions and period of analog simulation, 
permits the analog run to be completed, stores the re- 
quired information from the analog simulation, performs 
the required intermediate calculation, and proceeds on to 
the next reactor stage or tray. Problems of this type re- 
quire a prohibitive amount of time if only the digital is 
used and are too large for simulation if only the analog 
is used. 

Other hybrid applications have been primarily oriented 
towards the solution of mathematical relationships and 
application of mathematical techniques commonly en- 
countered in engineering. Solution of boundary value 
problems,'* application of the state variable concept"® and 
Monte Carlo method,*®*! statistical studies and 
correlation, *** and implementation of optimization 
schemes"®:*4#° are examples of the wide range of applica- 
tion and applicability of the hybrid in engineering. In all 
cases the use of the hybrid is highly recommended over 
either analog or digital computation alone. 

Of particular interest to the chemical and petroleum 
engineer is the use of the hybrid in solving partial differ- 
ential equations. Certain of these equations can be reduced 
to ordinary differential equations by an appropriate 
change of variable such as the familiar Blausius equation 
describing the laminar boundary layer on a flat plate® 
and numerous physical systems involving diffusion. 
Frequently a boundary value problem which is amena- 
ble to analog solution is obtained. In other cases the sepa- 
ration of variables technique which has direct analog 
utility may be used.* 

The solutions of other partial differential equations on 
the analog require finite differencing of all but one inde- 
pendent variable since only one independent variable can 
be continuous simulated on the analog. The use of this 
technique with only the analog can require an excessive 
amount of equipment although two independent variables 
can be simulated’. The alternative of digital simulation 
can lead to excessive computer time due to the large 
number of discrete values which must be handled in soly- 
ing these equations. 

"The hybrid affords the opportunity of time sharing the 
solution of partial differential equations resulting in a sig- 
nificant savings in analog equipment and time. The simu- 
lation of reactors previously noted’*"* represents one such 
use. The time required to simulate a two dimensional 
equation on a hybrid is substantially reduced over the time 
required for digital simulation since finite differencing 
need be used for only one independent variable. For third 
order equations, one independent variable can be repre- 
sented by analog operating time, a second by finite differ- 
ence approximations on the analog such as is done when 
solving second order equations on the analog,® and the 
third by finite difference calculations on the digital. 

The hybrid possesses the power which permits the 
solution of very large systems. Of importance to chemical 
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and petroleum engineering are applications such as model 
building for large systems,** the simulation of the wet end 
of a paper manufacturing machine, and the simulation 
of a chemical plant in order to train operators and opti- 
mize design.%°:33? 


The modern electronic analog computer is an in- 
herently powerful and rapid computational tool. The 
addition of analog memory and patchable digital logic 
capabilities establishes the analog as an important tool in 
engineering for the forseeable future. Despite the ad- 
vances in digital software, the marked speed advantage 
of the analog permits very large and complex engineering 
problems to be simulated in a reasonable amount of time. 

The combination of the analog and the digital, with 
appropriate interfacing, yields the most powerful and 
flexible computational tool available today, the hybrid 
computer. 
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How to Simulate Recirculation 


Reactor With Analog Computer 


Assuming recycle-fresh feed 
concentration for analog simulation 
proved reactor design was limited 


by turbulence level 


P. G. Cooper, Gulf Research & Development Co., 
Pittsburgh, Pa., and 

E. F. Schagrin, International Minerals & Chemical Corp., 
San Jose, Calif. 


For MANY HIGHLY exothermic reactions, it is necessary 
to dilute the feed to a commercial reactor with a recycle 
stream for temperature control, Often, means for deter- 
mining the initial composition to the reactor are nut 
available; thus, an initial condition problem is encoun- 
tered in the mathematical description of the process. In 
other words, the composition of the fresh feed is known, 
but the unknown composition of the recycle-fresh feed 
mixture is needed for the solution of the simulation equa- 
tions, The schematic diagram in Fig. 1 illustrates this type 
of reactor. 

After the startup of a Gulf processing unit incorporat- 
ing a recirculation reactor, it was found that an unusually 
low conversion was being obtained. It was suggested that 
the reactor design was unsound since the feed was quickly 
circulated to the exit point and could “short-circuit” 
through ‘he reactor. However, low reaction rate due to 
poor agitation could also give low conversion. It was 
Gecided to check the soundness of the reactor design by 
simulating it, assuming that the agitation was good. If 
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the computer results indicated high conversion, then the 
reactor design was sound and poor agitation was causing 
the trouble. 

The feed to the reactor is a two-phase liquid-liquid 
mixture and the heavy phase catalyzes the reactions in 
the light phase. Experimental measurements in a labora- 
tory continuous flow stirred tank reactor had shown that 
the rate constants k,, kz, and ks are a function of the tur- 
bulence level. The rate constants increased with increas- 
ing stirrer speed up to a certain level. 

From the chemistry of the system, a simplified kinetics 
model can be determined. There are primarily three re- 
actions that occur. AB, is the desired product. 


ky 

A+ B——> 4B (1) 
ks 

AB + B———>-AB, (2) 
ks 

A—————>C (3) 


The rate of heat transfer is sufficiently high so that the 
reactor is essentially isothermal and the reaction rate is 


tere 


PRODUCT 





COOLANT 


Fig. 1—Shows schematic diagram of jacketed recirculation 
reactor. 








+[s0 a8] +[50 a8], 
Fig. 2—Shows analog circuit diagram for the simulation of a recirculation reactor. 
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a function of only the point concentrations, A mathemati- 
cal model for a plug flow reactor can then be derived 
easily. Assuming no change in density due to reaction: 





—<o@ = ky (A)(B) + ks (A) 2 
a = ky (A)(B) + kz (AB)(B) (5) 
~ SAB) _ — by(A)(B) + a (AB)(B) (6) 
2 sae) = — kz (AB)(B) 2 
MIN (a) S 


The analog solution of equation (4) through (8) is 
straightforward, but we must first have the proper initial 
compositions to use in the simulation. We know that pure 
Aand B are mixed with the recycle stream before the inlet 
to the reactor; and the composition of the recycle stream 
will depend upon the inlet composition. Thus, we must 
resort to a trial and error procedure wherein: 


1. The initial composition is assumed. 


2. The reactor is simulated. 





3. A more accurate initial composition is calculated 
using the recycle composition indicated by the simu- 
lation. 


The procedure is repeated until there is no further change 
in the inlet composition. 


The analog circuit diagram for the solution of equa- 
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Fig. 3—Shows that part of the analog circuit to halt simula- 
tion when the outlet of the reactor was reached. 
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Fig. 4—Shows how the initial conditions converged to a fixed 
value for a given set of parameters. 


tions (4) through (8) is presented in Fig. 2. Q is the frac- 
tion of fresh feed in the flow to the reactor. The subscript 
f and i refer to the fresh feed and the reactor inlet, re- 
spectively. In order to facilitate the trial and error pro- 
cedure, a circuit was devised to automatically halt the 
simulation in the “hold” mode when the outlet of the re- 
actor was reached. This is presented in Fig. 3. When the 
simulation automatically stopped, the concentrations of 
all the components and the next set of initial conditions 
were read from a digital voltmeter. The new initial con- 
ditions were set on the computer and the simulation re- 
peated. The convergence of the initial conditions as a 
function of trial number is presented in Fig. 4. Approxi- 
mately 10 trials were necessary to obtain values that con- 
verged sufficiently for our purpose. 

The final prediction of the reactor profile (ie., conver- 
sion) compares well with the anticipated reactor per- 
formance showing that the reactor design is basically 
sound and that the level of turbulence is too low. To con- 
firm this conclusion, the reactor effluent was piped into a 
stirred autoclave which increased the conversion as ex- 
pected. Future reactor designs will include a much higher 
agitation level. 
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Fig. 1—The simulation program shows the regeneration policy which will maximize total production. 


Find Optimum Cat Regeneration Cycle 


An analog simulation technique 
determines the best catalyst 
regeneration policy. As various 
processing conditions change, the 
simulation program modifies 


the operating guides 


Michael D. Harter 
Electronics Associates, Inc. 
West Long Beach, N. J. 


Tue BEST CATALYST rejuvenation policy to be applied to 
a group of continuous reactors is determined as a func- 
tion of current operating conditions. An analog simula- 
tion technique is described which will show the best time 
at which to begin regeneration and upon which reactor. 
Return is measured as the total production over the 
planning interval, considering the productivity of each 
reactor being lost during its regeneration. 

The computational approach is to simulate, at relatively 
slow speed, the actual plant operation and catalyst activ- 
ity under the slowly-changing conditions forced by op- 
erator varied flow rate and severity. These current run- 
ning conditions are used as a moving base for the 
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high-speed investigation of alternative regeneration poli- 
cies, with an oscilloscope providing an operator guide in 
the form of a predictive display. 

The automatic determination of the best reactor to be 
regenerated and its corresponding regeneration time 
reduces the dimensionality of the total problem so that 
the remaining variables can be manipulated by the sim- 
ulation operator to meet planning needs subject to re- 
source constraints. 

The purpose of this article is to demonstrate a flexible 
and workable method to solve this particular recurring 
problem in a manner which takes advantage of plant op- 
erators experience. 

Since the main interest is in the control of solutions to 
the equations for the model and the investigation of best 
scheduling policies as influenced by different operating 
conditions, a very simple plant configuration and reactor 
description is adopted such as is shown in Fig. 1. More 
complex plants would not greatly alter the amount of 
logic required but of course would require additional 
analog equipment. 


EXAMPLE SITUATIONS 


Case 1. Assume a crude mechanism such as shown in 
Equation (1). This represents the formation of the prod- 


kk; 
A>B>C a) 


uct B from feed A, with degeneration to some unwanted 
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Fig. 2—Thruput and severity can be altered by operator. 
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Fig. 3—The No. 2 reactor is rejuvenated after 12 days. 
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Fig. 4—The peak production for each reactor occurs at dif- 
ferent times. 
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byproduct C. Assuming no B in the feed stream, the steady 
state solution to the differential equations: 


dA/dt = (q/v) (Aj, —A) —k,A (2) 

4B/dt = (q/v)(—B) +k,A—k,B (3) 
epeventre 

i 
LS 
(+8) ®) 
eS Ak,® = (Ain) (k,8) (5) 
(+R) ~ (FRONTED 
Giere 


@=0/q (6) 


It can be shown that the product concentration B is maxi- 
mized when: 





= (hy) (ky +k) = VOT 





In the particular case with ki = 1, Omar = Vika 


The reaction rates are assumed to be related to catalyst 
activity and severity of operation by the following: 


k,=sa (8) 
ky = ak, = asa (9) 

Instantaneous production rate is simply: 
’=4B (10) 
Thus, using Equations (5), (7) and (10), the best level 
of severity may be chosen once production rate is fixed 


or vice versa. 


Case 2. Now consider the following kinetic equation: 


k 
A>B (11) 
The steady-state solution to the equation 
dA/dt = q/v (Aj, — A) —k,A (12) 
is given by 


Differentiating the production rate of B with respect to 
flow and substituting for © gives 


d(qB)/dq = 4in( un y (14) 
q+ usa 





It can be seen that there is no supporting relation be- 
tween flow rate and severity which maximizes production 
in this case, and the variables (q and s) remain as inde- 
pendent degrees of freedom to the operator .For this rea- 
son the mechanism of Equation (11) was used. The in- 
stantaneous production rate as a function of flow and 
severity, is shown in Fig. 2. 


A word here about relative time-scales. The speed of 
response of a reactor to a change in operating condi- 
tions—either flow rate or severity—is of the order of 
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minutes. Catalyst decay time constants are of the order 
of days and hence the justification for the use of the 
phrase “instantaneous production rate.” To a first ap- 
proximation therefore we can treat the conversion kinetics 
as algebra when considered in relation to the planning 
interval of, say, one month. 


Degeneration. Decay is governed by severity of opera- 
tion and net production rate, the latter term correspond- 
ing to a poisonous influence. 

Tt is assumed that catalyst activity decay rate is 
increased by higher operating severities and is also propor- 
tional to the production rate. 


da/dt = —C [sa+ p) (15) 


This latter term corresponds to the deposition of carbon 
for example or gradual catalyst poisoning. 


Rejuvenation. It is assumed that rejuvenation takes place 
at a fixed rate from the activity level at which rejuvena- 
tion was begun up to a standard “clean” level. 


a —a (t,) 


T 


(16) 


r= 





Since activity decays in an exponential manner and 
rejuvenation normally takes place at low levels, the re- 
juvenation interval T, changes little if r is fixed. Thus 
the time for which production is lost is approximately in- 
versely proportional to 7. During the regeneration, begun 
at time t,, catalyst activity is raised from its level at that 
time up to some maximum level a? at the rate r. The time 
taken to regenerate 7’, depends on the time regeneration 
is started via a(t,), tending asymptotically to the maxi- 
mum a°/r as regeneration is postponed. 


PROCESS OBJECTIVE 


The approach to the maximization is a simple search 
over the space (i,j) for each value of a, q, and s. The 
variable « is predictive time extending from now to the 
end of the planning interval T. Note that the base r = 0 
is moving. The objective function then is: 


max (7) = 
ii) 


Ef Py (L(t) 9 (Fat eles (7-4 (Petty ae 
ate 

(17) 
with T<tyoct+T 
and OTA 


t, = j(T/32),j=0,1,2.- 





The accumulated product from each reactor over the 
planning interval for one particular regeneration policy is 
shown in Fig. 3. Note the plateau in the curve for reactor 
No, 2 and the improvement in production due to fresh- 
ened catalyst. These curves are for uniform operating con- 
ditions both throughout the interval T and between re- 
actors but for different initial catalyst activities 


(ay = @°, a, = 2a°/3, a, = 0°/3). 


As regeneration time is varied 72 (T) will rise and fall. 
If (T — T,-) <t, <T, the benefit from improved activ- 
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Fig. 5—Catalyst activity increases sharply after rejuvenation. 








Fig. 6—The simulation of three reactors in parallel. 
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Fig. 7—The analog circuit for reactor No. 2. 


ity is not felt, and an actual decrease in m2 (T) results 
from the shut-down. 


The influence of regeneration time on production is 
shown for each reactor in Fig. 4. The decreasing effect on 


FIND OPTIMUM CAT REGENERATION CYCLE . . . 














Fig. 8—Production from all three reactors is estimated to determine the best regeneration policy. 


the right is due to the shut-down loss just mentioned. 
Other points to be noticed are that the peaks for each 
reactor occur at different times and that A, the production 
improvement—the difference between = (t,) and x (7’)— 
is greatest for reactor No. 2. The rightmost points 7; (T) 
represent production if regeneration is not made within 
the planning interval. By inspection, the best policy in 
this case is P (2,4). 
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For policy P (2,12) and uniform operating conditions, 
Fig. 5 shows the three catalyst activities, initially a) = a°, 
a, = 2a°/3 and a, = a’/3. The slopes of decay at any time 
depend on operating conditions. For example, by reducing 
the severity on a reactor, it is possible to prolong the cata- 
lyst “life.” The area under the curves, a rough measure 
of production, is easily seen to be affected by a change in 
severity, throughput, and regeneration. For constant oper- 
ating conditions and a given regeneration policy, the pas- 
sage of time represents a continuous re-definition of the 
origin of o. 


COMPUTER IMPLEMENTATION 


The over-all organization of the simulation for a plant 
with three reactors in parallel, is shown in Fig. 6. 

This approach requires either a hybrid computer or an 
analog computer with sufficient parallel logic to imple- 
ment the search control, timing, and best policy selection. 
The solution requires many integration runs for the differ- 
ential equations of the models, needs convenient man- 
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Fig. 9—This method of search is used by the high-speed 
program. 


machine interfaces, and must run much faster than real 
time. Thus, an all-digital solution is ruled out at the pres- 
ent time. 

The plant is represented as three reactors in parallel, 
each running against the time variable 7. This is the low- 
speed section and accepts the operator input qj, si and t,;. 
This part of the simulation might also be used for opera- 
tor training. 

The high-speed section samples conditions in each re- 
actor in turn to determine the best policy at that time. 
This section runs against the time variable o, with 7 as a 
moving base. During the search, display of the individual 
7 as a functon of ¢, is available such as shown in Fig. 4. 
Following the search, the best policy may be displayed, 
showing either individual productions (Fig. 3) or indi- 
vidual activities (Fig. 5). 


The analog circuit for a single reactor is shown in Fig. 
7. Operator inputs are made through the potentiometers 
as indicated. The rate term to the “a” integrator switches 
between the decay variables—upper connection—and the 
regeneration variable r. Regeneration is initiated through 
a logic control push-button and is automatically termi- 
nated when activity is back to maximum. 

The high-speed circuit shown in Fig. 8 is multiplexed 
between the three reactor models just described. The only 
basic difference between this and the low-speed circuits 
is in the control of the regeneration switch. If the reactor 
currently being copied is in the regeneration mode, then 
all high-speed exploration runs for that reactor are ini- 
tiated in that mode also. This is in addition to regenera- 
tion as part of the search. 
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Fig. 10—Total production for a given policy is determined. 

































































Fig. 11—This circuit decides the best over-all policy. 


Control for the search mode is shown in Fig. 9. For 
three reactors and thirty-two days in the planning interval, 
a total of ninety-six runs are required for the exploration. 
Including initiallization times an exploration may be com- 
pleted in 200 milliseconds, an average 2.08 milliseconds 
per integration. This represents a scale-up of about 10° 
over real time. 

The circuit arrangement for computing zi, is shown 
in Fig. 10. The electronic switch in the input to the pro- 
duction integrator serves as high-speed hold during regen- 
eration. The lower right-hand three track/store units are 
used merely to display = (j) during search. The upper 
right three track/stores sample the production without 
any regeneration during the planning interval 7;. The left 
three track/stores pick the peak production for each 
reactor 77%. 

When the runs are completed, the circuit shown in Fig. 
11 decides the best over-all policy. The logic for this is as 
follows: 





If S > S, then Ly = 1. 
If S; > Sz then L, = 1. 
If $2 > So then L, = 1. 








FIND OPTIMUM CAT REGENERATION CYCLE . . . 


If Mo = (Lo) (L2) =1 
Tf Mi = (Ly) (Lo) =1 
If Mz = (Lz) (Li) =1 


then So best. 
then S; best. 
then S., best. 


The résults from decreasing the regeneration rate 25% 
are shown in Table 1. A comparison is made against the 


TABLE 1—Sensitivity to Regeneration Rate 














TABLE 4—Equipment Used 





ANALOG 
Amplifiers — 12 Track/Stores 
— 7 Integrators 
—17 High Gain 
— 7 Summers 
— 4 Electronic Switches 
— 38 Inverters 
85 
Multipliers — 12 (4 as Division) 
Pots —15 
Limiter ae 
Comparators —10 
D/A Switches — 22 
LOGIC 
And Gates —35 
Monostables —1 
Differentiators — 3 
Flip-Flops =$ 
Shift Registers — 2 
Mode Timer —1 
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base case of otherwise uniform operating conditions with 
an initial catalyst activities distribution of a = a’, a, = 
2a°/3 and a, = a’/3. As expected, the net effect is to 
decrease the maximum yield from each reactor. The best 
policy is still obvious although time of regeneration is ad- 
vanced. One useful application of this information might 
be to provide quantitative data for a regeneration de- 
bottleneck study. 

The corresponding data for a change to severity are 
shown in Table 2. The best policy under the new condi- 
tions (increased severity) is to regenerate reactor No. 2 
but postpone for two days the time of regeneration. On 
the other hand the recommendations with decreased se- 
verity is to advance the regeneration time two day. The 
change is production for this reactor is (0.4 + 0.1) and 
(0.4 — 0.12) respectively. 

Sensitivities to throughput changes are shown in Table 
3. The alteration in production for reactor No. 2—still the 
one to regenerate—is the same as in the case for severity 
changes, but in the reduced flow instance, an advancement 
of 4 days in regeneration time is necessary. 

Not shown in the table is the influence on 7;”: It is, 
respectively, 0.06; 0.04; and 0.02. This indicates that 
when reactor No. 2 is shut down, the available feed should 
be diverted unequally to reactors Nos. 0 and 1, with re- 
actor No. 0 having the larger share. Whether all the avail- 
able feed should be used and in what precise ratio it ought 
to be split must be determined at the time when reactor 
No. 2 is about to be shut down. 

All these sensitivity studies were made at a particular 
value of time r. By increasing the number of runs at each 
7 from 96 to 480, the sensitivity trials could be automated. 
The conditions fed to high-speed circuit would be: 


a(z),1,5(t), 4 (r) 

a(z), 1,5 (7), + As, q (r) 
a(z), 7,5 (1), — As, g(r) 
a(r),1,5(t),q(r) + Aq 
a(r),%5(t),q (7) — 4g 


For all i, j 


Time to run: 1 second. 


The corresponding changes to production and regenera- 
tion policy could then be displayed to the operator who 
may then decide to alter operating conditions as 7 evolves, 

Another alternative would be to use these influence 
coefficients as part of some higher-level automatic optimi- 
zation, say, steepest ascent manipulating s and q. 


NOMENCLATURE 


specific volume 
reaction rate ratio 


A,B,C component concentration v 
a 
A production improvement 
° 


a catalyst activity 


@° maximum catalyst activity 
k reaction rate constant residence time 
P regeneration policy total production 

: @ predictive time 
’ production rate FS 

7 present time 

q flow rate 
1 rejuvenation rate Subscripts 
S$ severity i reactor index 
T planning interval j day of rejuvenation 
t time 7 regeneration 


Indexing Terms: Analogs-10, Computers-10, Controlling-4, Optimization-4, 
Output-7, Performance-7, Processing-9, Production-7, Reactors-9, Regenera- 
tion-6, Simulation-8, Time-6. 
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This Distillation Program 


Generates Its Own Data 


This rigorous method incorporates 
the Chao-Seader correlation and the 
Redlich-Kwong equation of state to 
compute the equilibria and enthalpies 
for multicomponent calculations 


Wendell W. Waterman, Vern E. Alden Co., Chicago, 
and James P. Frazier, Natural Gas Pipeline Co. of 
America, Chicago 


Tue procra here presented is useful for hydrocarbon 
mixtures over a wide range of compositions, temperatures 
and pressures. It combines two important features: 


© Rigorous computation of the number of theoretical 
plates required for a specified separation. 


© Incorporation of the Chao-Seader correlation and 
the Redlich-Kwong equation of state enabling com- 
putation of reliable equilibrium and enthalpy values 
without introduction of data specific to each problem 
considered. 


PREVIOUS COMPUTER PROGRAMS 


Computer programs developed to solve multicom- 
ponent distillation problems by rigorous methods may be 
classed as follows: 


Class 1—Programs based on the Lewis-Matheson’* ap- 
proach in which feed, reflux and the desired separation 
are specified. The required number of theoretical plates 
is found by plate to plate calculation, component meshing 
at the feed plate being obtained by over-all iteration. 
Programs developed by McIntire,** Maddox and Er- 
bar'® and Bonner®° are examples of the Class I group. In 
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their use equilibrium and enthalpy data specific to each 
problem undertaken must be provided as part of the in- 
put data, Plate temperatures are found by computing 
bubble points below and dew points above the feed. 


Class M—Programs based on the Thicle-Geddes* ap- 
proach in which feed, reflux, quantity and state of the 
distillate, and the number of theoretical plates are speci- 
fied. Plate to plate computations are made using mol 
fraction ratios, and product compositions are found by a 
meshing computation. 

Programs reported by Lyster et al,* Mills,!? Newman"* 
and Hanson are of the Class II type. Again, enthalpy 
and equilibrium data must be supplied for each problem. 
However, in these programs determination of plate 
temperatures is more difficult. It is accomplished by an 
iterative process, a temperature profile being assumed 
to initiate computations. Various special techniques have 
been used to facilitate convergence. 


Class Ill—Programs based on solving the simultaneous 
equations which express plate and column material and 
enthalpy balances. The items specified are the same as 
for the Thiele-Geddes plate to plate method. 

The program developed by Amundson et al*:?8 belongs 
in Class III. Matrix inversion is used for solution of the 
simultaneous equations. Plate temperatures are deter- 
mined by iteration as in the Class IT programs. 


Class IV—Programs in which the progressive change in 
plate and product compositions from start-up to the 
attainment of steady state are computed by a relaxation 
method. Rose, Sweeny and Schrodt® have successfully 
developed this approach. 

All these programs are rigorous with respect to material 
and enthalpy balance computing techniques, but fail to 
take into account the composition dependency of equi- 
librium values. Further, although the tedious process of 






The method includes a rigorous computation of the number 


of theoretical plates needed 


making plate to plate calculations manually is eliminated, 
there remains the time consuming work of selecting and 
preparing equilibrium and enthalpy data for individual 
problems and the attendant possibility of human error. 

Greenstadt, et al’' described a program of the Class 
I type in which Newton’s method is utilized to solve 
plate and overall balances as well as the equations for 
adjusting product distribution to accomplish feed plate 
meshing. In one method of applying the program the 
Benedict, Webb, Rubin equation’ is used in a subroutine 
to compute equilibrium and enthalpy values. This was 


TABLE 1—Acceptable Components 








Non-Hydrocarbons Diolefins 
Hydrogen Propadiene 
Nitrogen 1, 2—Butadiene 


Carbon Dioxide 1, 3—Butadiene 
Carbon Monoxide 
Hydrogen Sulfide 


Sulfur Dioxide 


Cycloparaffins 
Cyclopentane 


Oxygen Methyl Cyclopentane 

Z Ethyl Cyclopentane 
Paraffins Cyclohexane 

Straight Chain—C, through C, Methyl Cyclohexane 

iso-Butane Ethyl Cyclohexane 

iso-Pentane ‘ 

neo-Pentane Aromatics 

2—Methyl Pentane Benzene 

3—Methyl Pentane Toluene 


2, 2—Dimethyl Butane Xylenes—Ortho, Meta and 
2; 3Dimethyl Butane Para 
ote Ethyl Benzene 


Straight Chain—Double Bond 
in 1 Position—C; through C; 
2—Butene—Cis and Trans 
2—Pentene—Cis and Trans 
2—Methyl—1—Butene 
2—Methyl—2—Butene 
3—Methyl—1—Butene 














TABLE 2—Input Specifications 








Feed 
Component List (Key Components Designated) 
Mole Fraction of Each Component in Total Feed 
Physical State Alternates 

(1) Liquid at Bubble Point 

(2) Vapor at Dew Point 

(3) Mixed Liquid and Vapor 
Temperature (For Flash Calculation) 


Distillate 


Physical State Alternates 
‘Vapor 
(2) Liquid 
Quantities of Light and Heavy Key Components 
Temperature (Starting Assumption) 


Bottoms Temperature (Starting Assumption) 


Reflux— 
Alternates: 
(1) Minimum Reflux Computed—Multiplier Specified 
(2) Reflux Ratio Specified 


Column Pressure 
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the first reported attempt to incorporate a generalized 
equation of state into a distillation program. It is under- 
stood that difficulty in attaining convergence has been 
experienced with this program. 

After an investigation of alternate methods for pre- 
dicting K values in hydrocarbon mixtures, Cavett’ 
adopted the Chao-Seader correlation® together with the 
Redlich-Kwong equation of state’® (or alternately for 
low pressure ranges the Black equation of state®) as the 
preferred method for obtaining thermodynamic proper- 
ties for use in distillation computations. Cavett developed 
two programs,® one of the Class II and the other of the 
Class III type in which it is necessary to specify the 
components in the feed, the reflux and the operating 
pressure. Equilibrium and enthalpy values are computed 
as required. 

There is an important difference in utility between the 
Class I and the Class II and III programs. The latter 
are analytic; that is, they compute performance in a 
column having a fixed number of theoretical plates. They 
are well adapted to studies on existing installations, but 
are inconvenient for design applications since the number 
of plates required to effect a given separation can be 
determined only by submitting a series of problems and 
interpolating among the results. The Class I programs 
are advantageous for design work since they compute 
directly the required number of plates. They also de- 
termine automatically the optimum point for feed intro- 
duction; whereas, with Class II and III programs this 
must be found by trial and error, 

Further the Class II and III programs usually have 
a maximum number of plates that can be considered. 
With Class I programs the number of theoretical plates 
is unlimited. 


CHARACTERISTICS OF THIS PROGRAM 


The present program is of the Class I type. It is used 
to compute the number of theoretical plates required to 
accomplish a given separation by a modified Lewis- 
Matheson technique. The Chao-Seader and Redlich- 
Kwong correlations are incorporated together with con- 
stants for the acceptable components listed in Table 1. 
Any 21 of the listed hydrocarbons plus any four of the 
non-hydrocarbons may be specified in a particular prob- 
lem. Thermodynamic data needed for rigorous plate to 
plate calculations are generated by the program for each 
vapor and liquid stream in the system. Equilibrium values 
are computed as functions of composition as well as tem- 
perature and pressure. 

Criteria and equations are included for checking com- 
puted feed plate compositions and making the adjust- 
ments necessary to secure non-key component meshing. 
The reflux ratio may be specified or the minimum reflux 
may be computed and used together with an assigned 
multiplier. 
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TABLE 3—Output Data 








Column Operating Pressure 

Feed Temperature 

Feed Vapor: Quantity, Composition and Enthalpy 
Feed Liquid: Quantity, Composition and Enthalpy 
Top Product: Quantity, Composition and Enthalpy 
Bottom Product: Quantity, Composition and Enthalpy 


Feed Plate Composition Calculated from Top for Each 
Meshing 
Iteration 


Feed Plate Composition Calculated from Bottom 
Minimum Reflux 

Minimum Reflux Multiplier (Specified) 

Reflux Ratio (Specified) 

Reflux Quantity and Enthalpy 

Reboiler Temperature and Vapor Enthalpy 
Reboiler Heat Load 

Temperature of Vapor to Condenser 

Condenser Outlet Temperature 

Condenser Heat Load 


Compositions, Temperatures and Enthalpies on Each Plate (If 
Desired) 











The program is written in FORTRAN and has been 
demonstrated on an IBM 7070-10K computer. It can 
be readily adapted to other computers that utilize FOR- 
TRAN compilers and have adequate storage capacity. 


Input specifications are listed in Table 2. It is neces- 
sary to specify only the feed, the required separation, 
the reflux and the column pressure. The feed may be 
liquid at the bubble point, any vapor-liquid mixture, or 
vapor at the dew point. Output data are listed in Table 
33 


At present the program is set up for one feed and 
two products. Modifications to accommodate multiple 
feeds and draw-offs are being developed. 


EQUATIONS USED IN PROGRAM 


x) ) 
%% 


log 49 = log ¥,() 4 ©; log ¥,0) (2) 


Chao-Seader Correlation. 





log ¥,0) = Ay +A,/T,+ 4,7, +472 +47 (3) 
+ (45+ 4,7, +4,72)P, 
+ (Ag+ 4gT,) P/t —log P, 


log »,0) = — 4.23893 + 8.65808 T, —1.22060/T, (4) 
— 3.15224 T,3 — 0.025 (P, — 0.6) 
_ 4 (5—8)) 
bees =—“9 503 RT ©) 






Redlich-Kwong Equation 





2° — 22 — zP (B,2P + B, — A,?) — A,2B,P? (6) 
* . 

Ag= Y) deve Be = Y 7sBrs (7) 
izt ist 


The compressibility, z, is used to calculate ¢; in the 
Chao-Seader correlation as follows: 


(z—1) By 

hk = > —B,P 

08 #5 = F503, log (z — B,P) (8) 
etl. Ass 


B, BP 
eee: =k) Jog ( 1.04 
Bes\ <A 2) os ( ae 


Enthalpies. Enthalpies of vapor mixtures are computed 
from ideal gas component enthalpies as follows: 











Ho= /)) »,H° (9) 
in 
A? P 
r= ne— re | sass log (10+ ) rt o-:| 
By, z 
(10) 
Liquid enthalpies are computed as follows: 
Aye= HY — 
= Ay Ag+ 7; (24,+34,7,) (11) 
Te : 
+ (Ag+ 24;T,) (Pp) + AoP/? 
1.22 
+ 6, (8.65808 + es — 9.45672 | 
—1.8y; (5 — 8,)? 
Hb = J) xy (12) 


in 
Material, Enthalpy and Component Balances. Over 
the entire column the balances are: 
F=F¥4FL=D4W (13) 
H,YFY 4 HEFL — q, + 9, = HyYDY + HeeD" + HW (14) 
(Note either DY = O or Dé = 0) 
x,FU + y,FY = y,DY + x,Do + xW (15) 


If balances are written around the portion of the tower 
below the line g-g in Figure 1, the following equations 
are obtained: 


Ly1=VatW (16) 
Aya! Ly + 9 = nV + yh W (17) 
Lg = IV te (18) 


Above the line g’-g’ the following balances may be 
written: 


Vy, = Dy + DY +D" (19) 
Hy? Vy — 9¢ = Hp Liga + HgYDY + Hrd” (20) 


IV mn = XiLma + DY + 2,DY (21) 








The accuracy of the results depends upon the accuracy of 


the Chao-Seader correlation 





Fig. 1—At present the program is set up for one feed and two 
products. 


Non-Key Component Meshing. The following equa- 
tion, developed by Edison, ° is used to compute new as- 
sumptions in iterations for meshing non-key components 
at the feed plate: 


Xfi — *fa,i 


45 =1.5)| Shot re, (22) 








This is the relationship employed by Bonner®® except 
for the arbitrary 1.5 multiplier which tends to speed 
convergence. 


Minimum Reflux. Minimum reflux is computed by the 
Bachelor method. Equations and computational tech- 
niques used in this operation are not presented here in 
detail but are available in Bachelor’s original paper. 
The special procedures for subcooled liquid and super- 
heated vapor feeds have not yet been incorporated in 
the program. 


COMPUTING TECHNIQUES 


The sequence of computations is shown in Figure 2, 
the desired separation being specified in terms of heavy 
and light key components. 
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The feed is first tested to see if a flash calculation is 
required. If so, it is computed to determine the split and 
the composition of the vapor and liquid portions. Next, 
the bottoms temperature is determined by bubble point 
calculation and the bottoms enthalpy is computed. The 
amount of reflux is computed and the overflow estab- 
lished at one point, usually the bottom plate in the 
column, This permits determination of the quantity of 
reboil vapor, its composition and enthalpy being available 
from the bubble point calculation. 

From this point on, material and enthalpy plate bal- 
ances are computed in the conventional manner, After 
each plate calculation the ratio of key components is 
checked to see if the feed plate has been reached. Com- 
putations of the condensing temperature, top product 
enthalpy, etc., are performed in a similar manner. 

For approximate calculations components heavier than 
the heavy key are assumed to be entirely in the bottoms 
and components lighter than the light key are assumed 
all in the distillate. The number of tlcorctical plates is 
computed in a single pass requiring from 7 to 12 minutes 
depending on the number of components and the number 
of plates. For most problems this procedure predicts the 
correct number of theoretical plates and determines ac- 
curately condenser and reboiler heat loads. 











If greater accuracy is desired the meshing feature of 
the program may be employed. The computer will adjust 
the amount of each non-key component in the top and 
bottom products until the feed plate composition as 
calculated from the top down agrees with that computed 
from the bottom up within a specified tolerance. The 
Edison’ algorithm used for predicting new assumptions 
as to distribution of non-key components has been found 
to converge in 4 to 7 iterations. Meshing satisfactory for 
most purposes is usually attained in four iterations. 

For each distillation problem a large number of bubble 
and dew point calculations must be made. These are 
trial and error computations, and to apply the Chao- 
Seader correlation it is necessary that compositions of 
both liquid and vapor phase be available. Bubble points 
are calculated by assuming a temperature, computing 
@ vapor, composition based on ideal K values, normalizing 
and applying the Chao-Seader equation to compute 
a new vapor composition. If SKjx; does not equal 
1.0 + the specified tolerance a new temperature is as- 
sumed and the process is repeated. After the second 
computation, when two sets of assumed and calculated 
temperatures are available, a linear fit is used to predict 
the next assumption. This simple algorithm has been 
found to converge rapidly and reliably. Similar con- 
vergence techniques are used for dew points, and for 
plate mass and enthalpy balances. 














The compressibility, z, is computed by trial and error 
using Newton’s method, the initial assumption being 
z= 1.0. In most problems difficulties due to extraneous 
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roots arise only at temperatures below the range of 
reliability for the Chao-Seader correlation. A rigorous 
analytic solution of the cubic equation has been pro- 





grammed, but is not normally used because of the in- a 
‘ cuannt 
Saran 


crease in computer time. 


Accuracy of Computations. The computation methods 
embodied in the program are completely rigorous and 
the precision of results depends only on the accuracy of 
the Chao-Seader correlation and on the tolerances set 
for convergence routines. Before the program was for- 
mulated Chao-Seader K values were compared with ex- agg ase 
perimental data and with values predicted by the Bene- varon aunty | 
dict, Webb, Rubin equation for mixtures of light paraffins 
and nitrogen. The following general observations were 
made: 























© The average deviation of computed K values from 
experimental values is about 8.0%. 


© The Benedict, Webb, Rubin equation predicts equi- 
librium ratios with slightly better accuracy at higher 
temperatures than does the Chao-Seader correlation. 


© At lower temperatures the Chao-Seader method 
is slightly more accurate than the Benedict, Webb, 
Rubin equation. 


The Benedict, Webb, Rubin equation, in numerous 
instances, will predict very accurately K values of 
one or two components of a mixture, but be con- 
siderably in error on other components. The Chao- 
Seader method gives more uniform accuracy with 
respect to all constituents. 








© The Chao-Seader correlation gives somewhat poorer 
predictions for binary and ternary mixtures than for 
mixtures containing a greater number of compon- 
ents. 


© Over its range of applicability deviations from ex- 
perimental data of K values computed by the Chao- 
Seader method are of the same order of magnitude 
as discrepancies among independent determinations 
by different experimenters. Fig, 2—The program is limited to ranges of conditions over 
which the Chao-Seader correlation is valid. 





Of the several K value correlation techniques that 
have been devised some are not easily adapted to com- 
puter applications. Some do not permit determination of 
equilibrium ratio as a function of composition as well as Experience. The program has been used to design dis- 
temperature and pressure. For general use none appear 
more accurate or versatile than the Chao-Seader cor- 
relation. © Natural gasoline stabilizer 


tillation columns of the following types: 


Program Limitations. The program is limited 10 ranges © Deethanizer 
of conditions over which the Chao-Seader correlation is 


Seria BITES aun oe Te Seo eee STE Coen pontents = eo heavier hydrocarbons 


and the number of components handled in a specific © Separation of iso and n butanes 
problem are described under an earlier section, “Char- 
acteristics of This Program.” © Separation of iso and n pentanes 


70 


DISTILLATION PROGRAM GENERATES OWN DATA 


© Debutanizer 
© Depentanizer. 


In each application many problems were computed 
both to test the program and to optimize column designs 
with respect to reflux ratios and operating pressures. 

During development, counters were provided in the 
program to record the number of iterations in each con- 
vergence routine. Various algorithms were tested and 
only those showing rapid and reliable convergence were 
retained. With the developed program no convergence 
failures have been experienced to date. 


NOMENCLATURE 
coefficients for simple fluids, methane or hydro- 
gen in Eq. (3) 
Ay, (0.4278 T,;25/P,,T?5] 05 
B,; 0.0867 T,,/P,;T 

D_ moles of top product 


0 








F moles of feed 
H molal enthalpy 

H® molal ideal gas enthalpy for mixture 
K, equilibrium ratio for component i 

moles of liquid 

P absolute pressure 

P, reduced pressure = P/P,,; 
P,, critical pressure of component i 

qe heat removed in condenser 

q, heat added in reboiler 

R_ universal gas constant 

T absolute temperature 
T,, critical temperature of component i 
T, reduced temperature T/T; 

V_ moles of vapor 
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W moles of bottom product 
mole fraction of component i in liquid 


mole fraction of component in feed plate liquor 
calculated from bottom up 


*fo,i 


ja, mole fraction of component in feed plate liquor 
calculated from top down 


y; mole fraction of component i in vapor 
z compressibility factor 

Y; activity coefficient of component i in liquid solu- 
tion 

adjustment to amount of component i in W for 
new assumption in meshing iterations 


solubility parameter of component i (in Chao- 
Seader equation) 





5 2x,¥,8,/2x,¥;, average solubility parameter of 
mixture 


»,° fugacity coefficient of pure liquid component i 
at system conditions 


¥,(0) fugacity coefficient of simple fluid in liquid state 
»,() fugacity coefficient correction factor 


$; fugacity coefficient of component i in vapor 
mixture 


¥; liquid molal volume of component i 
®; centric factor for component i 


Superscripts (other than those specifically defined) 
L liquid 
V vapor 


Subscripts (other than those specifically defined) 


d= top product 
f feed 
i any component 
m number of any plate above feed plate numbering 
from top down 
n number of any plate below feed plate numbering 
from top down 
w bottom product 
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Fig. 1—Key ratio in feed at optimum feed point matches key ratio in column. Sp = Sty = 





Compute Best Distillation Feed Point 


This computer program for 
multicomponent distillation is modified 
to give the optimum feed point 

for many problem situations 





Wendell W. Waterman, The Katzman Co., 
Youngstown, Ohio; James P. Frazier, Natural Gas 
Pipeline Co. of America, Chicago, and 

George Martin Brown, Northwestern University, 
Evanston, Il. 


Major ADDITIONS AND MODIFICATIONS have been made 
to a multicomponent distillation computer program* 
which broaden its range of applicability and enable it 
to cope with various special problems. The original pro- 
gram (now designated Mode I), based on a modified 
Lewis-Matheson? technique, computes the number of 
theoretical plates required for a specified separation, the 
feed composition being given and the operating reflux 
ratio being either specified or determined by application 
of an assigned multiplier to the minimum reflux ratio 
computed by the program. Data generation in situ is 
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accomplished by use of the Chao-Seader correlation and 
the Redlich-Kwong equation of state. 


The new features include: 


® A mode of computation (Mode II) that accepts 
problems in which the separation attainable with a fixed 
number of theoretical plates is computed, given the feed 
composition, reflux ratio, operating pressure, and the 
over-all split in terms of total moles of distillate and 
bottoms. 


© A mode of computation (Mode III) which accom- 
plishes preliminary determination of required theoretical 
plates by Mode I followed by meshing iterations under 
Mode II. 





COMPUTE BEST DISTILLATION FEED POINT... 


© A mode of computation (Mode IV) similar in essen- 
tial procedures to Mode I but which determines the 
optimum plate for introducing the feed in cases where 
the ratio of key components at the optimum entry point 
differs appreciably from the key ratio in the feed. 


© A mode of computation (Mode V) which accom- 
plishes preliminary determination of required theoretical 
plates based on optimum feed point location by Mode IV 
followed by meshing iterations under Mode II. 


e A routine for computing the tower operating pres- 
sure corresponding to a specified condensing temperature. 


© Amplification of the feed analysis routine to enable 
acceptance of superheated vapor and sub-cooled liquid 
feeds, 


© Storage of all constants used in the Chao-Seader 
and Redlich-Kwong equations for the entire list of 60 
acceptable components in a block data sub-program. 


© Major simplification of problem entry techniques. 


Mode Il. In the previous article,’ it was explained that 
the Edison modification’ of the Bonner equation’? was 
used for computing new assumptions with respect to 
distribution of non-key components in successive meshing 
iterations. Experience has shown that the Edison modifi- 
cation, involving an arbitrary multiplier applied to the 
computed adjustment for each component, tends to speed 
convergence but in certain instances causes excessive 
swings in early iterations. Consequently, the Bonner 
equation is now used in its original form. 


X jw, — Xa, 


Xya,i 
TDi 





To initiate a problem, it is the usual practice to enter 
as the starting assumption for the mole fraction of each 
light component in the bottoms and each heavy com- 
ponent in the distillate, a very small figure, 1 xX 10-*°. 
Not infrequently, this amount proves too great, and as 
the computation proceeds from plate to plate, excessive 
increase in the mole fraction of one or more components 
may occur. 

This “overbuilding” distorts the computation and makes 
it fruitless to continue without altering the assumed dis- 
tribution. Occasionally, too, a value predicted by the 
Bonner equation for either a light or heavy component 
will be too great and will lead to overbuilding. This is 
particularly apt to occur in the second meshing iteration. 


The program incorporates, both above and below the 
feed, check routines which compare at each plate the 
computed mole fraction for each component with the 
match figure it should agree with (or cross over) at the 
feed plate. If the match figure is greatly exceeded (some 
tolerance is allowed), and if the key ratio indicates that 
more plates must be computed before reaching the feed 
point, new assumptions are set for distribution of the 
overbuilding components and computation of the tower 
section under consideration is restarted. Further, once 
a value for a component in the distillate (or bottoms) 


REPRINTED FROM HYDROCARBON PROCESSING 










FEED COMPOSITION s 





An 


OPTIMUM: 
FEED POINT 








BOTTOM - 


PLATES TOP. 


ratio at imum feed point is below that 
of teed. Soon = Eaaaty 3, oe 


is found to be excessive, this is established (with a toler- 
ance to allow for composition dependency interactions) 
as a limit so that subsequently no appreciably higher 
value will be tried even if called for by the meshing 
formula. Similarly, lower limits are established when 
values too low for meshing at the feed plate are tried. 

The Bonner equation guided by the overbuilding 
checks and the progressively narrowing limits proves to 
be a powerful tool in handling meshing problems. It is 
found that problems involving a fixed number of plates 
above and below the feed may be fed directly to this 
routine and solved reliably with rapid convergence. It is 
necessary, in submitting a problem of this type, to assume 
an initial split that represents the desired quantities of 
distillate and bottoms respectively, 

This is the method employed in Mode II. The Thiele- 
Geddes approach* has been widely used for solution ol 
problems of this type in other computer programs.'-* The 
present technique offers the advantage that it does not 
require initial assumptions as to either temperature or 
composition profiles over the column. In Mode II, key 
components are not distinguished in the computations 
and all components are simultaneously meshed at the 
feed plate. 


Mode Ill. This Mode was developed to cope with an 
unexpected difficulty encountered in certain problems. 
In successive meshing iterations under Mode I, the dis- 
tribution of non-key components in the top and bottom 
products are adjusted progressively until a pre-determined 
tolerance is met or until a limiting number of meshing 
iterations has been run. In most cases the number of 
theoretical plates is not affected by changes in com- 
ponent distribution from one iteration to the next. How- 
ever, in certain instances, an oscillation in the number 
of plates either above or below the feed, though usually 
not both in a single problem, may occur in the following 
manner: 
Assume that m plates above the feed have been found 
for the Nth meshing iteration based on a set of assump- 
tions as to amounts of non-key components in, for ex- 
ample, the distillate. For the N + Ith iteration new 
assumptions are computed. Results of the N + Ith 
iteration show m + 1 plates required to reach the 
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Fig. 3—Key ratio at optimum feed point is above that 
of feed. Sb, 4 St,. ee 





proper key ratio match point. Further, the proportions 
of heavy components on the feed plate are found to 
be excessive since the amounts in the distillate were 
set for meshing on the mth plate and computations 
were continued past this point. The N + 2 assumptions 
are then computed on the premise that meshing must 
be accomplished on the m + Ith plate. This calls for 
reduced amounts of heavy non-key components in the 
distillate, the amounts being similar to those assumed 
in the Nth iteration. Upon running the N + 2 itera- 
tion m plates above the feed are again found and the 
computed quantities of heavy components on the feed 
plate are too low. Oscillation in the number of plates 
will continue until an assigned limiting number of 
meshing iterations has been reached. 
It might appear that this oscillation could be avoided by 
slightly altering the specified split. However, in cases 
encountered so far, this was not possible. Increased 
stringency in the separation merely shifted the oscillation 
point, so that, for example, instead of oscillating between 
m and m + 1 plates, the computation oscillates between 
m + 1 and m + 2 plates. 

Integer plate number oscillation occurring in this 
fashion is due to relative variation in the K values of 
key components stemming from composition changes. It 
conclusively demonstrates the necessity in truly rigorous 
computations of using composition dependent equilibrium 
values. 

Difficulties arising from plate number oscillation might 
be avoided by determining the fraction of a theoretical 
plate required for key ratio matching at the feed point. 
However, another method of coping with the problem 
has been devised and incorporated in Mode III. This 
Mode may be used if oscillation is discovered in a prob- 
lem submitted under Mode I. 

In Mode III two iterations are run under Mode I 
and the larger number of plates found, above, and also 
below, the feed is accepted and fixed, and operation is 
transferred to Mode II. This leads, in some cases, to a 
slightly more complete separation than was ori inally 
specified with respect to key components but provides a 
practicable and rigorous solution to the problem. 

Mode III can, of course, be used directly for problems 
in which no oscillation is either discovered or suspected. 
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Fig. 4—At low reflux ratio an extremely large number 
of plates below feed is required to reach feed key ratio. 


It has the advantage, since key designations are dropped 
after two iterations, of meshing all components at the 
feed plate and thus more accurately represents actual 
tower performance than does the Mode I computation. 


Mode IV. Problems in which the optimum feed point is 
displaced from the plate on which the ratio of key com- 
ponents matches that of the feed have been encountered 
with sufficient frequency to warrant development of a 
feed point optimizing technique. 

For binary systems, the compositions of vapor and 
liquid leaving the optimum feed plate are directly related 
to the feed composition. As a first assumption, in multi- 
component systems, the ratio of key components in the 
feed vapor (or liquid) may be taken as equal to the key 
component ratio in the vapor (or liquid) leaving the feed 
plate. This criterion gives reasonably good results in 
many cases. It has been observed to represent roughly 
the optimum feed point when the relative volatility be- 
tween the key components is significantly closer to unity 
than either the relative volatility between the light key 
and the next lighter component or the relative volatility 
between the heavy key and the next heavier component. 
When only one or neither of these conditions exists, the 
key ratio on the optimum feed plate frequently has been 
found to be appreciably displaced from the ratio of keys 
in the feed. Other factors may also affect the location 
of the optimum feed point. 

Plates computed from the top down and from the 
bottom up are plotted in Fig. 1 on independent ordinates 
against the key component ratio in the liquid for a 
typical multicomponent problem in which the feed, the 
tower pressure and an operating reflux ratio well above 
the minimum are specified. In the case illustrated, the 
optimum feed point falls at the point where the key 
ratio in the vapor from the plate matches that in the 
feed vapor (or the vapor portion of the feed, or the 
vapor in equilibrium with a liquid feed). It will be noted 
that in computing plates from the top down the key 
ratio falls rapidly at first, then at a diminishing rate. 
Finally when carried far beyond the key ratio of the 
feed, the key ratio asymptotically approaches a limit rep- 
resenting a very large number of plates, and has an ex- 
tremely small change from plate to plate. 


COMPUTE BEST DISTILLATION FEED POINT... 


A similar, but inverted curve results in the computa- 
tion of plates from the bottom up, a large change in 
key ratio is experienced at the beginning and an ex- 
tremely small change as an asymptotic limit is ap- 
proached. The spread between the asymptotic limits 
represents the range of possible solutions to a specific 
problem. With hydrocarbon mixtures such a range typi- 
cally exists for any problem in which the operating reflux 
ratio is well above the minimum, The optimum solution 
is that which gives the minimum number of total theo- 
retical plates. 

When the key ratio on the optimum feed plate matches 
that of the feed, as illustrated in Fig. 1, the tangent lines 
to the curves at this key ratio will be parallel. It is 
apparent that the minimum number of plates for any 
particular problem will be found at the key ratio where 
the slopes of the tangent lines to the curves are equal, 
whether or not this key ratio coincides with that of the 
feed. 





These slopes are designated S'opr and Sopp in Figs. 1, 
2 and 3. Fig. 2 illustrates the case where the key ratio 
at the optimum feed point is appreciably below that of 
the feed. Fig. 3 shows the optimum feed point falling 
at a key ratio higher than that of the feed. It will be 
noted that in Figs. 2 and 3, the slopes of the tangent lines 
at the key ratio of the feed Sy and S®» are unequal and 
that the sum of top and bottom theoretical plates is 
greater than at the optimum feed point. In Mode IV, 
the program computes as much of each plate vs. key ratio 
curve as may be required to find the point where the 
slopes are equal and the number of theoretical plates is 
minimum. 

The use of this optimizing technique is facilitated by 
entering, for the first iteration, small numbers for the 
assumed mole fractions of heavy components in the dis- 
tillate and light components in the bottoms since this 
permits optimization for the first iteration to proceed 
with minimal chance of overbuilding. However, if needed, 
the overbuilding routines will operate to enable the 
computation. Location of the optimum feed plate, as 
found in the first iteration, is subject to redetermination 
in subsequent iterations, but such redeterminations are 
not apt to cause a shift of more than one theoretical 
plate. After each computation the key ratio on the opti- 
mum feed plate is stored and used as the match point 
for the next meshing iteration. 

At reflux ratios close to but above the minimum a 
number of variations in shape of the plate vs. key ratio 
curves may occur. At the true minimum reflux ratio only 
one point of feed introduction gives a solution and the 
ratio of keys on the feed plate may not coincide with 
the feed key ratio. Consequently it is possible to en- 
counter problems which, over certain ranges of reflux 
ratio, have real solutions but which present no solution 
based on feeding the column on the plate where the key 
ratio matches that of the feed. 


A problem is illustrated in Fig. 4 in which a normal 
curve exists for top plates but the curve for bottom 
plates approaches the feed key ratio very slowly. In this 
case an extremely large number of plates below the feed 
is indicated when closing on the key ratio of the feed. 
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With the optimizing technique the best solution, repre- 
senting minimum plates for the separation specified and 
the reflux ratio considered, is readily found, 

An example is presented in Fig. 5 in which a hook 
exists in the top plate curve making it impossible to reach 
the feed composition key ratio. The bottom plate curve, 
however, continues well beyond the key ratio of the feed. 
In executing the program this condition is detected by a 
slope inversion test. When inversion is found top plate 
computations are stopped, bottom plate computations are 
initiated and the point of parallel slopes is found in the 
usual manner, A problem of this type cannot be solved 
by direct convergence on the key ratio of the feed but a 
range of solutions nevertheless exists and the optimum is 
determinable by the program. Fig. 6 illustrates the same 
type of problem but with the hook occurring in the 
bottom plate curve at a key ratio below that of the feed. 


A low reflux ratio problem is illustrated in Fig. 7 
which exhibits another variation in the interrelations 
among plates, meshing and possible feed points. When 
the problem is introduced in the usual manner with very 
small amounts assumed for heavies in the distillate and 
lights in the bottoms the first meshing iteration gives an 
indicated optimum of the type shown in Fig. 2. How- 
ever, with new assumptions respecting the distribution 
of non-key components, in the second iteration the re- 
sults show hooks in both the top and bottom plate curves 
with no overlap of key ratios. 

In this particular problem, any attempt to distribute 
non keys to secure meshing at the key ratio where this 
first set of equal slopes is encountered will result in dou- 
ble hooks with no common key ratio. Despite indications 
at this point, possible solutions may exist, as shown, the 
entire range falling between the inversion points of the 
top and bottom curves respectively. 

The optimum feed point is found at the point where 
the slopes of the inverted curves become equal. The zone 
of possible solutions may be appreciably less than the 
full span between low and high invert points since it is 
limited to points at which complete meshing of all com- 
ponents is possible. 

Problems of this type are rarely encountered. When 
they are, practical considerations usually call for higher 
reflux ratios. Consequently no attempt has been made 
to include in the program the capability of optimizing 
and meshing in this special situation. If such a problem 
is submitted for solution, the existence of non-overlap- 
ping double hooks is detected and the computation ter- 
minated. The type of problem presented in Fig. 7 can 
be resolved, however, by use of the program, Mode II 
will compute the meshing of all components in this case 
with key ratios following the curves as shown. However, 
the number of total plates required and the optimum 
feed point for a specified separation and given reflux 
ratio must be determined by trial and error in successive 
entries of the problem. 

Typical curves of plates vs. key ratio at minimum 
reflux are illustrated in Fig, 8. In this case only a single 
mesh point is possible at one particular key ratio, which, 
as indicated above, may be either above or below the key 
ratio of the feed. 


Computation of Column Pressure. In the original 
program, column pressure was specified and the con- 
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densing, reboiling and all plate temperatures were com- 
puted. In many cases it is desired to specify condensing 
temperature in order to meet a given cooling water or 
refrigeration temperature level. The new routine for 
computing tower pressure entails a determination of the 
condensing temperature at two assumed pressures and 
calculation of the tower pressure by interpolation or 
extrapolation based on an assumed linear relationship 
between the logarithm of the absolute pressure and the 
reciprocal of the absolute temperature. 

If the condensing temperature is computed on the 
basis of an initially assumed distribution of components 
in the distillate, a slight variance may occur when a more 
precise determination of component distribution is made 
in meshing iterations. Usually the change is not enough 
to warrant re-computation of the operating pressure. 


Acceptable Feed Conditions. In the original version 
the program accepted only saturated vapor, mixed liquid 
and vapor or saturated vapor feeds. If the specified feed 
temperature was found to correspond to a sub-cooled 
liquid or a superheated vapor, the specified temperature 
was ignored and the feed was taken at its saturation 
temperature. The new version of the program accepts 
any feed temperature and determines whether this cor- 
responds to liquid, vapor, or a mixture of liquid and 
vapor. If a superheated vapor or a sub-cooled liquid 
condition is found, the enthalpy is computed and the 
feed is introduced at the specified temperature. Enthal- 
pies of liquid, vapor and liquid-vapor mixtures are com- 
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puted by equations developed by Edmister, et al."* The 
vapor enthalpy equations are applicable to vapor in 
both the saturated and superheated conditions. The 
liquid enthalpy equations were developed for the bubble 
point liquid and have not yet been tested for sub-cooled 
liquids. In the program these have been used, however, 
for sub-cooled liquids; so for feeds in this state computed 
enthalpies must, at this time, be regarded as approxi- 
mate. 


Block Data Storage of Constants. It has been found 
convenient to store all constants required for the Chao- 
Seader and Redlich-Kwong equations for the entire list 
of acceptable components in a block data sub-program. 

th the original program, it was necessary to submit 
with each problem a previously prepared card, for each 
component in the feed, on which appropriate constants 
were punched. Block data storage eliminates this require- 
ment. 


Problem Entry. Problems can now be entered by pre- 
paring one card for each component plus two additional 
cards. On each component card the component number 
and name, and identifying number are punched together 
with mole fractions in the feed, distillate and bottoms 
and an appropriate designation if the component is a 
light or heavy key. In addition to the component cards, 
one card is required for problem identification and an- 
other for selection among alternate computing modes. 


In Modes I, III and IV, with the components listed 
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in order of decreasing volatility, heavy components in the 
distillate and light components in the bottoms may be 
specified as zero or left blank on the problem input cards. 
If meshing is called for, the program assigns to such 
components the starting value 1 X 10.*° The split of 
key components must be specified by entering an amount 
for the heavy key in the distillate, and for the light key in 
the bottoms. If for either key the other amount is omitted 
the program computes it. A starting assumption with 
respect to split must be provided for distributed com- 
ponents falling between the keys. 


The input formats permit entry of any desired set of 
starting assumptions with respect to distribution of non- 
key components in the distillate and bottoms. Thus, if 
reliable approximations are available in advance they 
may be used to minimize the number of meshing itera- 
tions. Further, if a problem has been run for say three 
meshing iterations, the computed distributions may be 
entered as input in re-submitting the problem, in which 
case, meshing iterations are continued without repeti- 
tion, 

In the case of hydrogen in the bottoms the initial 
assumption of 1 X 10-*° has been found too high in 
many instances. Any desired value may be entered and 
this will override the arbitrary assignment. It is the usual 
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practice to enter 1 X 10-*° for hydrogen in the bottoms, 
particularly if a fairly large number of plates below the 
feed is anticipated. 

A convenient problem input form has been devised 
consisting of a single 8% X 11-inch sheet on which all 
necessary information may be entered. The form includes 
all required instructions to the key punch operator so 
that the time of technical personnel required to set up 
and enter a problem is minimal. Once defined, a prob- 
lem can ordinarily be entered in 15 to 30 minutes. 


Computer Time. Core time has been considerably re- 
duced by numerous minor improvements since the de- 
velopment of the original program. Over the last year, 
many problems have been run on IBM 7090/94 and 
Control Data 3400 computers. In this application the 
speeds of these machines are about equal. A 17-compo- 
nent Hydrotreater Effluent Stabilizer problem with com- 
ponents ranging from hydrogen to decane has been 
found to require about 21 minutes when run through 
five meshing iterations. Satisfactory meshing of all com- 
ponents was achieved. In the IBM 7094 computer the 
core requirement including that of the processor pro- 
gram and library subroutines is approximately 18 K, 36- 
bit words. 


NOMENCLATURE 
D_ moles of top product 
A; adjustment to amount of component i in W for new 
assumption in meshing iterations 
K, equilibrium ratio for component i, (y;/x;) 
ky moles of light key component in liquid 
ky moles of heavy key component in liquid 
m number of theoretical plate above feed counting from 
top down 
N_ number of meshing iterations 
Sb, slope of key ratio versus plates below feed curve at 
key ratio of feed liquid 
Sty slope of key ratio versus plates above feed curve at 
key ratio of feed liquid 
slope of key ratio versus plates below feed curve at 
optimum feed point 
slope of key ratio versus plates above feed curve at 
optimum feed point 
W moles of bottom product 
*, mole fraction at component i in liquid 
ya, mole fraction of component i in feed plate liquor 
calculated from top down 
mole fraction of component i in feed plate liquor 
calculated from bottom up 
9, mole fraction of component i in vapor 
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Analog Simulation Spells Safe Startups 


Analog simulators have long been used 
in operator training for normal 
operating conditions. But training for 
vital startup phase has been largely 
ignored. Here is how the Humble-Esso 
organization bridged the gap 


R. E. Lieber, Esso Research & Engineering Co., Florham 
Park, N.J., and T. R. Herndon, Humble Oil & Refining 
Co., Baytown, Texas 


In AN EFFoRT to provide experience before startup, 
the techniques of dynamic, realistic simulation have been 
used in a Esso Research & Engineering Co. training pro- 
gram for experienced operators. These techniques involve 
simulating operations of a highly integrated section of a 
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new refinery on an analog computer. The goal is to dem- 
onstrate major process interactions and their effect on 
plant operations. 

Plant type instruments tied into the computer serve as 
the “control house.” This training program has proved 
highly successful, as this refinery has been brought on- 
stream smoothly. Thus, more training sessions of this kind 
are planned. 

Startup simulation is especially useful to an operator 
on a new, complex, highly heat integrated plant. He has 
a unique problem because his plant has never before been 
started up and is inherently more difficult to operate 
than older plants. 

There is benefit in extending such a simulation to cover 
startup—that is, bringing the plant up to design operat- 
ing conditions. This would be of value not only in terms 
of improved operator training, but also in setting opti- 
mum startup procedures. 

There are big differences between a model developed 
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FIGURE 1—This section of the refinery was simulated by the 
computer and model (Six distillation towers, with preheat 
exchangers and reboilers, and two reactors). Notes: Unit con- 


for normal plant operation and one that simulates start- 
up operations. After describing the normal operation 
model, we will point out these differences. We will con- 
sider what is involved in a plant startup and how this 
affects the computer model. 


Analog Model Adds Realism. The section of the re- 
finery simulated for operator training is shown in Figure 
1. Included are six distillation towers, with preheat ex- 
changers and reboilers, and two reactors. For over-all heat 
economy, 11 exchangers transfer heat from one process 
stream to another, rather than to cooling water. The 
plant has no intermediate tankage. The only storage con- 
sists of four 15-minute holdup drums with the outflow of 
each on flow control. 

A mathematical model was developed for each major 
piece of process equipment, describing its steady-state 
and transient behavior. The complete plant was pro- 
gramed on an analog computer containing about 400 am- 
plifiers and 400 potentiometers. This plant simulation, 
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LEGEND 
R = Reactor 
T = Distillation Column 
C = Primary Temperature 
Controllers (controllers 
and recorders on 
control console) 
Indicates Unit Heat 
Integration 








tains no intermediate tankage. Not shown: Four 15-minute 
holdups (flow control on outflow); All other associated equip- 
ment (condensers, coolers, etc.) ; Pressure, flow, level controllers. 





occupying four standard consoles of analog computing 
equipment, was then connected to a control console. 

The control console has standard plant recorders and 
electronic controllers. Their input and output voltages 
are compatible with the computer. 

The analog model and control console together simu- 
lated the operator’s normal plant environment. It also 
gave him access to similar records and indications that he 
would have in an actual control house. To use time most 
efficiently in the 10-day program, simulation was time- 
scaled, speeding up the process by a factor of 10. 

Recorder speeds were stepped up by the same factor. 
Thus, recorder charts appeared familiar to operators. By 
taking advantage of the analog computer's time-scaling 
feature, many dynamic control situations could be inves- 
tigated in a short time. 

Trainees were all experienced operators, So emphasis 
was on familiarizing them with operating and control 
characteristics peculiar to the unit. The model gave di- 
rectional indications of system behavior under upset con- 
ditions. It also illustrated the time relationships involved. 
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Thus, operators got the feel of the instruments. They 
learned how various disturbances propagate through the 
plant, what effects their corrective actions have, and how 
in a complex situation these actions feed back through 
the processing units. 

Training sessions pinpointed control problem areas 
(such as where automatic control is essential) and con- 
ditions when operator intervention is required. 


Time-Scaling Compresses Experience. An introduc- 
tory period was devoted to familiarizing operators with 
equipment in use. The period included a short descrip- 
tion of the computer. The time-scaling factor built into 
the simulation received special attention. It was illus- 
trated with several simple control loops. 

Instructors then demonstrated the use of the model. 
They showed how to identify upset conditions, how to 
anticipate required control actions, and how the heat in- 
tegration circuit is affected by disturbances. 

At this point, trainees took over operation. Typical dis- 
turbances were introduced into the simulated plant. 
Trainees studied the effects of these disturbances and the 
effects of corrective actions. After the simulations, results 
were analyzed in group discussions. Typical disturbances 
considered included changes in feed rates, compositions, 
set-points, feed preheat temperatures and reactor outlet 
temperatures, and major upsets such as pump failure and 
complete loss of instrument air. 

The differences between manual and automatic con- 
trol were illustrated. These are important in a highly in- 
tegrated unit. The trainees learned what the automatic 
controllers can do without operator assistance. They were 
taught what an operator can do to improve control under 
abnormal conditions. 

In the final training phase, the instructors applied dis- 
turbances which, on the basis of their experience with 
the model, the operators identified and corrected for. 


Simulating Process Units. Even with the large amount 
of computing equipment used, certain simplifications in 
the mathematical model were required. Since operator 
training was the primary purpose, operability limits on 
some of the equipment were not considered. 

For example, it was assumed that the fractionating 
towers can hydraulically handle the design rates or any 
applied deviations. Therefore, the model did not include 
dumping and flooding correlations. Since pressure con- 
trol is relatively fast and usually good, constant pressure 
at the design value was assumed for the model. 

Control console instruments were mainly recorder-con- 
trollers for temperature and differential temperature at 
points in the fractionator towers, pre-heaters, and else- 
where in the heat integration circuit. Such local instru- 
ments as level indicator-controllers on low holdup vessels 
were not simulated. Other recorders (i.e., flow) were sim- 
ulated by 8-channel analog recorders. 

In deciding how to simulate the dynamic behavior of 
refinery equipment, the purpose of simulation is very im- 
portant. For this program, the main concern was to build 
into the simulation directional indications that familiar- 
ize the operators with the idiosyncracies of particular 
units, Errors up to + 10 percent were tolerated. 


80 


Tem. . ts) ke tS 
Borie Varo Vis) Tse ays 


























k 1 2 6 
T a e v 
ave J O fs 
7 
r e | 
" 7 60. } 


FIGURE 2—Analog model strikes a balance between realism. 
and computer capacity, Depending on type of column and 
its complexity, model has 7-20 amplifiers and 10-25 poten- 
tiometers. Distillation column simulation is shown here. Trans- 
fer function relating control tray temperature to reboiler vapor 
obtained from digital model. 
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FIGURE 3—Computer model for kettle reboiler simulation 
(constant shell side temperature). 


The model was to represent daily operations. So, ex- 
treme excursions from design points would be rare. This 
allowed certain simplifying assumptions, such as linear- 
ized perturbation models, tower control temperatures as 
functions of composition only and not of pressure, and 
pseudobinary tower feeds. 

To simulate major items of refinery equipment we em- 
ployed existing rigorous digital models, from which 
linearized transfer functions, were derived, These trans- 
fer functions were then approximated with simpler ones. 
The simplifying assumptions were based on previous op- 
erating experience with similar units. In some cases, ac- 
tual plant data from a similar unit was used as a check. 

Because heat integration affected the reboilers and con- 
densers, distillation columns were isolated from the over- 
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head and bottom equipment in the calculations. Dynamic 
plate-to-plate calculations were made digitally, assuming 
pressure to be constant (but not necessarily uniform 
across the tower). 

The feed was broken down into pseudobinary compo- 
nents, 

Linearized transfer functions were obtained to relate 
product compositions, overhead and bottoms tempera- 
tures, and control tray temperatures to such inputs as 
preheat, reboiler heat, reflux and feed rates, and product 
withdrawal rates. In simplifying these transfer functions, 
a balance was struck between program realism and com- 
puter capacity. Depending on the type of column and 
its complexity, the analog model, Figure 2, contained 7 
to 20 amplifiers and 10 to 25 potentiometers. 


Realistic Reboiler Simulation. Simple heat exchangers 
such as kettle reboilers with constant shellside tempera- 
ture and no phase change in the tubes were simulated 
rigorously. They required 7 amplifiers and 10 potenti- 
ometers. For the more complex exchangers, the dynamics 
were first derived digitally, then approximated by sim- 
plified transfer functions. 

Each heat exchanger was broken down by tube sec- 
tions. For each section, heat and material balances were 
programed to solve for tubeside inlet temperature to the 
next section @, and shellside vapor generated V. 

The analog model for a kettle reboiler, Figure 3, solves 
the equations: 

do, 
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where the subscripts T, S, and M refer to tubeside 
fluid, shellside fluid, and tube metal; h is heat 
transfer coefficient, A heat transfer area, © tempera- 
ature, m mass, and c specific heat. A is the latent 
heat of vaporization of the shellside fluid (tower 
bottoms) ; wr is the mass flow rate of the tubeside 
fluid (reactor effluent), and my is its holdup mass 
in the ith section of the exchanger. 





Reactor simulation required extensive steady-state cal- 
culations. Corr-lations were first developed relating out- 
let temperature and conversion to changes in inlet tem- 
perature, flow, and concentration, Based on flow rate, 
holdup, and assumed flow patterns, dynamic behavior 
was calculated. It was then approximated by such simple 
transfer functions as a dead time and two time constants. 
Reactor simulation circuits were similar to those for dis- 
tillation columns. They required 3 to 6 amplifiers and 
4 to 6 potentiometers. 

Temperature transmitters were simulated by a first-or- 
der lag. The dynamics of transmission lines, valves, 10 to 
15-second holdup vessels, and flow controllers were ig- 
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nored. These responses can be assumed instantaneous 
relative to the important control aspects of the unit- 
temperature control and liquid level in 15-minute holdup 
drums. 


Startup Operation. Shortly after the training program 
was completed, the refinery was brought onstream. Al- 
though some of the usual mechanical problems occurred, 
the plant was smoothly lined out. No problems arose in 
areas the simulation covered. 

Similar analog simulations for operator training are 
now underway, not only for highly integrated plants, 
but also for new processes where no previous operating 
experience exists. 

As plants become more complex with higher heat and 
process integration, the startup phase becomes more criti- 
cal. With the success of the program just described, we 
investigated how to simulate actual startup. There are 
several distinct phases after construction of a new plant 
is completed. These involve the pre-startup activities such 
as commissioning utilities, running in of pumps and com- 
pressors, flushing lines, drying furnace refractory and 
steaming out towers. 

After the unit is ready to go, circulation is started by 
pumping streams into various pieces of equipment. 

Levels in drums are established, Heat is added to warm 
up the equipment and process streams. Inert gases are 
vented and pressures are built up, The startup progresses 
as additional equipment is brought into operation in 
order and the unit begins making product. These streams 
are initially recycled or sent to slop. The final stage of 
startup consists of lining out the operating units to bring 
the products on specification. 


Startup Involves Crises. This description of a startup 
makes it appear a simple operation. Actually, it is a se- 
ries of small crises. A number of process and mechanical 
difficulties usually develop. These must be corrected be- 
fore the next step can be taken. 

Thus, the operator must know how to cope with 
equipment failures occuring at any time in startup. He 
must be ready for emergency conditions, know how to 
recognize them, and know what action to take. He 
must manipulate a large number of process variables 
either manually or by controller set point adjustments 
as the unit is brought up. Instruments are normally 
commissioned as the startup progresses. Still, the operator 
cannot heavily rely on his controls to correct upset and 
emergency conditions. 

In a complex heat integrated plant, the operators’ 
difficulties are compounded. With heat integration, the 
startup sequence will usually not follow the processing 
sequence. This is because the heat source for upstream 
equipment may not be available until downstream equip- 
ment is operable. Thus, pieces of equipment may have 
to start up simultaneously and be brought up together 
to maintain heat balance in the unit. 

Also, without intermediate tankage between units, 
levels in surge vessels occupy much operator attention. 
Thus, the startup of an integrated plant places a heavy 
load on the operator. This burden is lightened if he 








has been through the startup a number of times and is 
aware of what to expect from action he takes. 


The Same Computer Models Do Not Apply. A start- 
up simulation requires a mathematical model of the 
time varying behavior of each piece of equipment. 
However, the two models are vastly different. 

Simulation of normal operations involves small ex- 
cursions about an operating point. Certain variables 
such as tower pressures and levels in low hold-up vessels 
can be considered as constant. With this approach, 
wide use can be made of linearized perturbation models. 

During the startup of a plant, flows, temperatures, 
and pressures undergo many-fold changes. Few variables 
are constant for a computer model. Since startup is 
sequential, its simulation undergoes a major change as 
each piece of equipment is commissioned. Therefore, 
the main emphasis of a computer model is to account 
for the heat balance and flow of material as the start- 
up sequence progresses. 

Compositions of product streams do not become sig- 
nificant until the line-out phase. For example, in a 
pipestill, the initial startup step will involve gas oil 
circulation. This brings up temperatures in the tower 
and preheat exchanger train before crude oil is brought 
into the unit. With the changes in the process variables 
involved during startup, properties such as gas densities, 
heat capacities, and thermal conductivities also cannot 
be assumed as constant. They become dependent vari- 
ables in a computer model. Thus, perturbation models 
have little application and many nonlinear computing 
components are required. 


Simulating a Tower. As an example of the differences 
between the two models, consider a fractionating tower. 
For normal operation, about the design value, a very 
adequate model consists of establishing steady-state re- 
lationships and using linearized transfer functions to 
relate tower temperatures and product rates and com- 
positions to input variables. These transfer functions 
result in simple analog circuitry. 

For the startup of a fractionating tower, a computer 
model must consider warm-up, pressuring and purging. 
During the pre-startup phase, the tower will have been 
steamed out and pressured with inerts or light gas. As 
feed and reboil heat are introduced, hydrocarbons begin 
condensing on the tower surfaces. This is the principal 
means of heat transfer during the warm-up. However, 
the heat transfer coefficients depend strongly on the 
amount of noncondensible gas present in the tower and 
this changes as the tower is purged. 

Describing the transients of pressure, temperature and 
the condensed hydrocarbon accumulating in the tower 
and distillate drum results in nonlinear differential equa- 
tions such as the following ones for tower pressure. 
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where Pr, Py, and P; refer to totals, hydrocarbon 
and inert pressures in the tower; Ty and Ty are vapor 
and metal temperatures; V is tower volume, R is the 
gas constant, Qe is reboiler heat, 


) 
Py 
is the heat transfer coefficient as a function of the 
hydrocarbon composition, \ is latent heat of vaporiza- 
tion and O is the overhead vapor rate. 

Similarly, differential equations with dependent co- 
efficients are involved for the vapor and metal tem- 
peratures. The analog circuitry for the set of equations 


describing this operation is considerably more complex 
than that shown in Figure 2 for a single tower model. 


Can the Two Models Be Combined? Since the 
simulation of normal operation and the startup each 
require much analog computational equipment, the 
question of equipment limitations is raised. Obviously, 
a combination of the two models would be the most 
effective training tool. Mathematically, this is not dif- 
ficult. Simulating the operation and startup of a small 
and single unit is easily done. 

However, the value of computer simulation is in its 
application to complex and integrated processing units 
from which small portions cannot be isolated. Thus, 
the nature of the application demands extensive com- 
putational facilities. The problem that now remains is 
the development of solution techniques to permit a prac- 
tical meshing of the two models. 


This paper originally presented at the 1964 annual 
meeting of the Instrument Society of America, October 
12-15, New York City. 
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Train Power Station Operators by Analog 


American Oil’s Whiting refinery has 
used analog simulators successfully 
in a variety of training programs to 
upgrade operators. Here is their 
latest innovation, which you may use 
to increase employe proficiency 


L. G. Whitesell and E. H. Bowles 
American Oil Co., Whiting, Ind. 


A RELATIVELY NEW, interesting, fast and efficient 
method for using an analog model to train refinery 
power station operators is enjoying success at American 
Oil Co.’s Whiting refinery. Models are proving useful in 
many training applications as well as for simulations 
from which to make operating calculations. (For example 
see “Train Your Operators by Computer,” HP/PR, July 
1963, Page 192; and “Try This Analog for Solving Fire- 
Water Problems,” September 1962, page 274). 

The latest innovation is a computer-operated working 
model for power stations. It is set up to fit unusual oper- 
ating problems at Whiting refinery, but could be applied 
to any refinery where an analog computer is available. 
Here’s how the system works. 


A Complex System to Learn, The two separate, but 
tightly interconnected, steam and electric power gen- 
erating stations at the Whiting refinery make up a highly 
complex system. As electric demand has increased and 
process steam demand has decreased in recent years, 
operating problems have steadily grown more difficult. 

The efficient training of new men and upgrading of 
older men’s skills has become more important as the 
steam-electric power balance has grown tighter. 

Operators must be able to comprehend the overall 
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FIGURE 1—These two stations provide all of Whiting refin- 


ery’s electric power and most process steam. Line widths show 
relative magnitude of steam flows in various parts of system. 


system, so they can function smoothly as a “team” to 
achieve a common operating goal. 


The Whiting Power Stations. All electric power and 
most of the process steam required by the refinery is 
produced by two stations (see Figure 1). The No. 1 
Station produces 400-psi steam which is passed through 
turbines and converted to 100-psi steam for process use. 
The No. 3 Station produces 1,250-psi steam, which is 
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FIGURE 2—The system is controlled by this analog computer, built 
in a lab near the refinery. It is used exclusively for simulation’ work. 


passed through three turbines to produce 400-psi steam, 
and through one turbine to produce 100-psi steam. Three 
small turbines in the No. 1 Station (only one is shown 
in the figure) use 100-psi steam and provide added oper- 
ating flexibility in special situations. 

Waste heat recovery types of steam generating units 
are located on refinery process units and beyond the con- 
trol of the power station operators. These provide half 
the total refinery demand for 400-psi process steam. 

Pressure reducing and desuperheating units parallel 
each of the two main sets of turbines. They are auto- 
matically controlled to maintain a set pressure in the 
steam system downstream of the reducer. The desuper- 
heaters normally absorb the steam flow variations caused 
by changes in steam and electric demands. 

The line widths in Figure 1 show the relative mag- 
nitude of steam flows in the various parts of the system. 


Operating Problems. The station operators face two 
basic classes of problems. One class involves unusual 
demand changes and equipment failures, when operating 
conditions are tight and flexibility is limited. 

The other involves ability to minimize production of 
“non-byproduct power” during day-by-day normal opera- 
tion. 

In the past, process steam demand well exceeded that 
produced by turbines in meeting the electric demand. 
And, enough steam flowed through the desuperheater 
units to accommodate any sudden change in demand for 
either steam or electricity. Individual refinery processing 
units were small enough that sudden large changes in 
system load were unlikely. 

Today, steam demand barely exceeds that produced 
by the turbines alone, and the refinery consists of a 
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few large processing units. So, the likelihood of sizeable 
swings in demand is much greater. 


Byproduct Power. Power produced by the turbines 
when a process requires their exhaust steam is called 
“byproduct power.” In this case the cost of the fuel 
represented by the latent heat in the steam is charged to 
the process. 

Such byproduct power requires about 3,600 Btu per 
kilowatt-hour. But if condensing or atmospheric exhaust 
turbine drives are substituted for the usual electric 
drives on process units, the “non-byproduct power” thus 
produced could require as much as 45,000 Btu for every 
kilowatt-hour, 

One can see that minimizing the production of non- 
byproduct power will save money. Clearly, when the 
steam-power demand ratio is tight, the reward for opti- 
mum day-to-day operation is great, During such a period, 
much skill is required to handle some types of sudden 
demand changes, because minimum steam flow is desir- 
able through the desuperheaters. 


Development of a Trainer. Several years ago, an 
intensive operator training program was undertaken by 
American Oil to improve power station operating tech- 
niques. It included both classroom and on-the-job instruc- 
tion. The results were mostly satisfactory. But new opera- 
tors required something more than the usual forms of 
instruction because they lacked the years of experience 
of older operators. Simulation seemed a good answer. 


Computer Plus Control Board Equals Training. An 
analog computer facility designed to simulate refinery 
process equipment is located in a laboratory near the 


85 








refinery. Associated with it is a large control board 
equipped with commercial recording and control instru- 
ments. These can easily be arranged in a graphic-type 
display of the process simulated on the computer, which 
js monitored and controlled by the board instruments 
themselves. This installation, used to train operators for 
refinery process units, is also used to train power station 
operators. 

The power-station trainer features a graphic instru- 
ment display which follows, as closely as possible, the 
physical layout of equipment in the two stations. Record- 
ers and indicators are provided for all important steam 
flows, steam pressures, and electric power demands. Fig- 
ure 2 shows the electronic analog computer used to simu- 
late the system. 


Control Board Setup. Figure 3 demonstrates how a 
boiler, a turbine-generator, and a desuperheater are 
represented on the control board. The trainee has control 
knobs for starting up, shutting down, and trimming the 
outputs of individual boilers. Switches adjacent to these 
knobs permit the operator to switch from automatic to 
manual control of boiler firing. He can also choose low 
or high speed draft fan operation. 

Wherever possible, the simulator has been tailored to 
produce system responses which “feel” real to the experi- 
enced operator. The knobs controlling turbine governor 
adjustment, for changing the relative loadings on the 
turbine-generators, present the same action and “feel” as 
their counterparts in the actual stations. 

Even system frequency has been made to respond 
realistically to load changes and machine adjustments. 
Behind the board, compressed air is vented noisily when- 
ever header pressures exceed the settings on pressure 
relief valves. 

A small panel out of sight of the trainee is used 
by the instructor to vary the electric loads on feeder 


SSS 
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FIGURE 3—Here is how a boiler, turbine generator and de- 

superheater are represented on the simulator’s control board. 
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groups, change the demands for steam, and cause the 
failure of selected station and refinery equipment. 


The Training Program. Training on the analog simu- 
lator is conducted in two phases. Both are under the 
direction of an experienced shift supervisor. 

In the first phase, one boiler room operator and one 
turbine room operator are together familiarized with the 
simulator panel. Then they are given a standard set of 
typical operating problems to solve. Situations they must 
handle include 1. Emergency loss of a boiler under 
several conditions. 2. Loss of a turbo-generator. 3, Un- 
expected shutdown of a major refinery process unit. 
4, Minimizing of non-byproduct power generation during 
a period of tight steampower balance. 

This session lasts eight hours, The second phase is 
much like the first, However, the operators bring their 
own problems, questions, and experiences to the training 
session. All current trainees have completed phase one, 
and several have completed phase two. 


Results and Conclusions. Because the program is still 
under way, the full results of this training technique have 
not yet been evaluated. However, there are many gratify- 
ing indications that the program is working out well. 

‘Analog simulator training does not eliminate the need 
for classroom and on-the-job instruction, nor does it 
detract from the value of many years of experience. But 
it does provide a quicker, clearer means of getting a 
complex operating problem across to an inexperienced 
operator. 

The performance of the trainees and the questions 
they ask reveal the weaknesses and gaps in their under- 
standing and previous training that were not apparent 
before. These weaknesses point up areas in the con- 
ventional training program which need to be reinforced. 

Finally, turbine room and boiler room operators are 
becoming more conscious of one another’s problems. It 
is apparent that a better informed, more thoroughly 
rounded-out crew is being developed. ## 








Minimum Reflux Figured an Easier Way 


David L. Ripps 
American Cyanamid Co., Wayne, N.J. 


To caLcULATE minimum reflux ratio by the Under- 
wood method, it is necessary to solve the following equa- 
tion’ for 6: 


I 


Ou X+ 
2 Sees ~ 
U (®) y 285 (1+ p)=0 (1) 
t=1 
where a; = relative volatility of component ¢ 
x; = mole fraction of component # in the feed 
p thermal condition of the feed 
ZI = number of components 





Equation 1 is highly non-linear with a singularity at 
each value of © = aj. The roots of interest lie in the open 
interval between adjacent pairs of «’s. Because of the 
strongly non-linear character of Equation 1 a direct ap- 
plication of the standard methods of iteration (viz., New- 
ton, false-position, etc.) often fails to converge or requires 
an excessive number of cycles to reach a root. We will 
now develop a specific method for solving Equation 1 in 
a reasonably efficient manner. 


In this discussion it will be convenient to number the 
components so that the e’s are in ascending order. 


am <a <<... Say (2) 
Let us suppose that we wish to find a root, say the one 
in the interval 
a, <O < art 


Presuming for the moment that the root lies close to a, 
we may rewrite Equation 1 as 








Oy Xr 
i = ¢, (0) (3) 
where 
Oy X; 
¢, (8) = ean U(®) (4) 


Equation 3 may be solved for © by treating $, as a con- 
stant. Since ¢, is actually a function of ©, the solution 
must be iterative. Then 





(a = (5) 
where 
(Q), = guessed value of © at cycle n of the iteration 


(@.), = computed value of © at cycle n of the iteration 
Expanding Equation 5 in first-order Taylor series gives 
40. 
(@)a+1 — On = Sy [. s17 | 
A ssolution on the next cycle would require 
(@d)n+1 = ()n41 


Combining the foregoing two equations and solving for 
the next cycle 


(8.)n — (O)n 
1 — (d0,/d0) 

For efficient solution we will use the previous two 
cycles in order to evaluate the derivatives, 





(O)n41 = On + (6) 


dO. _ (8c)n — (Oc)n-1 
do (8)n — (O)n- 
To start the procedure we will take © as the mid-point 
between a, and a;.1. 








=2,3,... (7) 





(@): = 0.5 (a, + a, +1) (8) 


We assumed originally that the root in the interval be- 
tween a, and a;.; is closer to a, than to a,.;. Evaluating 
the function at (@); provides a test of this assumption. 
If U(®); is positive the assumption is correct and we 
may proceed. If that function is negative, however, the 
root is closer to ay,1 and we should replace r by r+1 in 
subsequent parts of the iteration. 

To continue we use Equation 5 as a direct iteration to 
give the second starting point. Thus 


(0): = (@.) 


Thereafter we will use the modified Newton given by 
Equations 6 and 7. 


Example Problems. 
Example 1: Benzene — Toluene — Xylene? 





m = 0.1 0.3 x = 0.6 
am = 04 = 1.0 as = 2.45 
p= =I 


Desired: Root in interval 1.0 < © < 2.45 


Cycle | i um 








1 | 1.725 1.5836 1.1502 
2 1.15020 — 0.9198 6 1.2784 
3 :25501 6 } =11834 1.2535 
4 1.25380 2.8228 x 10-3 | —1:1821 | 1,258 





Example 2. Four component system. 





Desired: Root in interval 1.0 < @ < 1.3 
Tolerance on 6 + 0.00001 


Cycles to 








* Race |. @ r_| converge 
1 005 | u10 | O75 | 0.10 | 1.02705 | 2 4 
2 0.05 | O10 | 020 | 0.65 | 1.05060 | 2 4 
3 025 | 025 | 025 | 0.25 | 111967 | 2 4 
4 o10 | o40 | o40 | 0. vig6er | 2 4 
5 065 | 0.20 | 0.10 rigo71 | 3 4 
8 a10 | 075 | v.10 | 0 125504 | 3 4 
ae eS) SSS 


LITERATURE CITED 
1 Robinson, C. S., and Gilliland, E. R., “Elements of Fractional Distilla- 
tion,” McGraw-Hill Book Co., Inc., 1950, p. 356. 
?Ibid., p. 358. 


Indexing Terms: Computations-4,7, Design-t, Distillation-8,9, Engineering-4, 
Equations-10, Iteration-6, Retluxing-9, Separation-8,9, Underwood-0. 








TABLE 1—Input Data 


































































































1 5 r 5 ape ea RL 
iia00 Ne oa = 1 ene po 22 5 ae 
if iia 1 + aoe ee oo oe SS ee 
Sot pa goat ee a ee a ee ES ad 3 es 
tt ot tt Sie ey See ee ts 
3 treo i a0 a iia ipa a st ae se ge 
ie oor sr et 3 ee 
Sina emer ee cs tietens Tyce se gp iay o oeae ee sre meee oie EET SY So ages 38 3 to 
pues ee ee ea ae et 7 toe od 
3 i000 ¢ ont ete 251 ay slaaiiee ao vinse. 0 Seana 3? wo 
Be ee = yeu —— 
i: Rage AN Et serge ek gm We ging Crag ea gy as tT 0.0 
+H tt Sa a te = oe iF a 
7 "gne,08 1 0.0, a ee ae 3 a 
$e et aot = 3 T a 
Sa a pnt Une eNO pos pay eo gaara pet 3 oo 
$B ae Oe ae te t — I 
3 e009 4 0.0 EME ST i Ser sis a 2.0 
St te r 3 iT 
Re iec Wane ere ge ge cg pe ag ese 3 ee 
Sercae tt oa $i et a OO é eee 7 
Rip eagaaire Dantes er ghee eRea  guear gy Se yeas e840 ea 
+ 4+ i a a : fe ot 
Sie ra Npeise sateen Siping gece ge ye ema pepe g Ab saestie 
Sars os eH : att t + — Seon 
Yai eee mas ni cokes WEEa DNs pags at ae ove een sy 1 wee 
ft Sa ee a ee 
Top eae nc ats as UROL Seg se se ys =e fei: ae 
tart i Bt oe 
ey ere mega eae eng see easy 4 iieeseaEeEis cg hf 2 4eee pp mgreen ry ta 
Sa ee ed a—t Bt tet 7 
SS ST a a ial cone reeey ze 903) 93315 oat 
hott et st + st 3 SI 
GT NSS ST es eS Pap aeat apreien os 
$i pH Ht i ot So 
og epg epee gs eb ag eases 3 
ehh tore qe 1 rt So 
Pee ee ares ee) igs eg Faye sess -e = 
jt et tt Pt Oot SLL Ea 
ES a ee a a 1 
tet I set ot i Sa a 
Later er ee ener ye Sag Ee Teg tome ib toe i 
#—t $d et : a a ar 
1 150.0010 0.0 ois oy priate toss : 
4 cr saa 3 aoer 
Pea eee Vee Pay paat. seear bs oe 0 6 8 0 
it ir tt | 1% 
1b 357.0048 0 oie 3 ool ke gS bi 321022428 
a se =x r aes + it 
Fea a 3 3 2: Fe gecenreene ip eued 
+ tt 5 a a ee 
2 ¢ a. a cone OE ae a gaa 
*5 3 Pome eye a et uo ee 
st ot 














Design Best Cooling Water System 


This computer program helps design 
an economical cooling water system. It 
gives the overall installed cost of the 
piping and can be used for sizing 
the pump 


Jogindra P. Kohli* 
Union Carbide Corp., Tonawanda, N.Y. 


Tuts COMPUTER PROGRAM is based on the installed cost 
of the piping for a cooling water system at a given loca- 
tion. For given allowable pump head, pipe routing and 
maximum and minimum velocity limits; it designs the 
most economical piping system. The program can also be 
used for sizing the pump, calculating head losses in an 
existing system, and computing the cost of a given system. 
Typical results from the computer program are presented. 


Optimum Design. What is a most economical design? 
It is the design which satisfies all the required conditions 
and costs the least. 

For a cooling water system, the required conditions to 
be satisfied are (1) flow, (2) maximum velocity, (3) 
available head. Flow through a pipe is equal to its area 
times the maximum velocity and therefore two of the 
three requirements can be satisfied by calculating the min- 


* Present address: State University of New York at Buffalo, Amherst, N.Y. 
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imum diameter of pipe for required flow and maximum 
velocity, The system thus designed is the best design if it 
can satisfy the third requirement of available head, but 
usually it does not. Therefore, some lines in the system are 
to be increased in size (within the range of minimum 
velocity) such that available head loss conditions are satis- 
fied. 

Which lines should be increased? The answer can make 
the difference between the best and the worst design. If 
the lines that are increased cost the least, the design is 
the most economical design. It is not an easy job to find 
which lines cost the least to increase (especially in a large 
system). For example, a line of a relatively smaller diam- 
eter with two valves may cost more to increase in size 
than a larger diameter line of some length but without 
the valves. This is where the computer can be used to the 
greatest advantage. 

The design of a cooling water system (or a water sup- 
ply system) involves a trial and error procedure. In the 
normal procedure, the lines are sized approximately for 
the given flow and the friction head loss is calculated. 
This friction loss is added to the various other head losses 
(equipment head loss, cooling tower head, etc.) and 
checked against the available pump head. If the head loss 
calculated for the system is more than available pump 
head, some line sizes in the system are increased and fric- 
tion head loss recalculated. The procedure is repeated 
until the total head loss through the system becomes less 
than the available head. 





TABLE 2—Results of Cooling Water Piping Design Program 





CIRCUIT Noe 
3 


HEADLOSS(FT 7 
740 

















10 
1 





i aL 
ALLOWABLE PUMP HEAD = 108.0 COOLING TOWER HEAD = 20.00 





A OST OF THE WHOLE SYSTEM IN DOLLARS = _29076,03 





DESIGN DATA USED_IN THE PROGRAM 
MAXIMUM ALLOWA IGN VELOCITY FLOM 45 GPM AND MORE = 12,00 ET./SEC. 
MAXIMUM ALLOWABLE DESIGN VELOCITY FLOM LESS THAN 45 GPH = 1U-00 FT. /SEC> 
AVEKAGE LABOR RATE FOR INSTALLEO COST CALCULATIONS = 7,00 S/Miinoun 


INSTALLED COST UNDER GROUND PIPING IN DOLLARS = 21366.98 
INSTALLED COST ABOVE GROUND PIPING IN DOLLAR 7709-96 


INSTALLED COST OF THE WHOLE SYSTEM IN DOLLARS = 24076103 a 





TABLE 3—Head Loss Through the Circuits 











BRANCH FLOW OIA WEADLOSS INSTALLER COST —AUIWE GATUND Uni a 
HO, G.PAM, EET poliaKs co COST WNL LANs) 








5557.00 14.000 Tess ‘s206-33 














Bele pule uf | 










8 
1-000 
42000 





Et 


21 630.50 6000 
22__05,00__¢:000__ 
23 2000 1.300 
-24__230,0¢ 100 
253921450 14.000 

28 


27 —125R;50 10-000 
fk 3 4000 
29 —~T93ns00 12,000 
30__952:00 102000. 
31 825.00 10-000 


3 05:00 __6:000. 
: 0 

34__3421:50__16:000, 
3am om 

36 2583100 __ 14:00. 





Types of Problems. For a given value of allowable 
pump head, a defined routing, given flow and limiting 
velocities the program designs the cooling water system 
for least cost and calculates the total installed cost of the 
system. The program can design a number of systems for 
different head values for the same data and thus can be 
used in selecting a pump, such that the combination of 
pump cost (installed cost plus operating cost for certain 
number of years) and the system cost is least. The pro- 
gram can handle a system up to 30 circuits (each circuit 
having not more than 48 lines) and up to 120 lines (each 
line having not more than 16 fittings). The program can 
calculate the head losses in an already sized system, and 
it can also compute the cost of such a system. 


Method of Analysis. The procedure consists of four 
basic steps: 

© Initial sizing of the system 

© Determining the head losses in the circuits and com- 
paring them with the available head 

© Determining the most economical line to be increased 
in size and repeating operations one to three until head 
losses in all the circuits are less than the available head 
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© Computing the cost of the system thus designed and 
printing out of the results. 

Initial Sizing of the System. Initially all the lines are 
sized to satisfy flow requirements with maximum velocity. 


Area = Flow/ Velocity 
xD'/4 = (144 (0.1337)Q)/(60 AMV) 
where 





D 
Q = flow through line, gpm 
AMV = allowable maximum velocity, ft./sec. 
Then 
D = 12 0.1337 Q/47.1 AMV a) 
The lines thus sized are the minimum sizes for given 
maximum velocity conditions. 
Head Loss Through the Circuits. After lines are sized, 
head loss through each line is calculated by means of the 


following two equations: 
For length of pipe (Hazen Williams Formula), 


(0.002083) (QZ) (QF)**/(PDc)""* (2) 


pipe diameter, in. 


AL 
Where 
HL = friction head loss, ft. 
QL = length of pipe, ft. 
QF = flow through pipe, gpm 
PDC = pipe diameter, in. 


W 


For fittings in the pipe, 

HF = (0.00259) (K) (QN) (QF)*/(PDC)' (3) 
Where 

HF = friction head loss, ft. 

K = constant for different fittings 

QN = number of fittings 


The head losses through circuits are calculated by add- 
ing equipment head losses to the line head losses. These 
circuit head losses are compared with the available head 
(pump head-static head). The circuits, with head loss 
greater than available head, are used to determine the 
most economical line to be increased in size to reduce 
the head. 

The assumptions are: 

© Labor rate is assumed to be a uniform average rate 
for a particular location. 

© Value of “C-100” is used for old piping at 60°F in 
calculating head loss in piping by William and Hazen 
formula. 

© Values of “K” for fittings, in calculating head loss 
by KV*/2g, are used as per “Hydraulics” by Schoder and 
Dawson. 


Most Economical Line. Pipe sizes are increased one at 
a time, and each time a most economical line is chosen. 
Selection of the most economical line is done in the fol- 
lowing way: 

Costs of all the lines that occur in the circuits requir- 
ing head loss reduction are calculated. All these lines are 
increased in diameter to their next size and their new 
costs are calculated. The difference in the two costs for 
each line gives us the amount of money it will cost to in- 
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crease that particular line by one size. The same opera- 
tion is repeated for head loss calculations and difference 
in two head losses gives us the amount of head loss re- 
duction due to increase in that line size. This is the actual 
head loss reduction. 

Effective head loss reduction is obtained by multiply- 
ing actual head loss reduction by a factor that depends 
upon the effectiveness of the increase of size of the line 
to the whole system. The factor known as the effective- 
ness factor is obtained for each line. To calculate the 
effectiveness factor for the line, the number of times 
each line occurs in the circuits requiring head loss is 
found, If the difference between the circuit head loss and 
available head is more than the head loss reduction due 
to increase in size of line, the value of factor is one for 
that circuit. Otherwise it is the ratio of the difference 
of circuit head loss and available head to the head loss 
reduction due to increase in line size. If the line does not 
occur in a circuit requiring head loss, its effectiveness 
factor is zero. All the factors for all the circuits are 
summed for each line. This gives the value of the effec- 
tiveness factor for each line. 
cost per foot of head for each line is obtained 
ing the amount it costs to increase that line one 
size by the effective head reduction due to that increase 
in size. These costs are compared for all the lines and 
only the line with the least cost per foot of head is 
chosen. That line is checked for minimum velocity re- 
quirements. If it satisfies, this is the most economical line 
size to increase. If not, the line with the next higher cost 
is chosen. 

The above chosen. line is increased in size and head 
losses recalculated. The procedure is repeated until all 
the circuits satisfy the head loss requirements. This is the 
optimum design. 

When all design requirements are satisfied, the pro- 
gram calculates the cost of all the lines and prints out 
the line sizes, head loss, cost of each line, head loss 
through each circuit, and the cost of the over-all system. 





Computer Program. The computer program has been 
written for IBM/360 computer which has a storage ca- 
pacity of 64,000 locations in the core storage. The pro- 
gram consists of a main program and five subroutines and 
requires 49,282 locations. Input data is read from the 
cards and the output is on a paper printer. The program 
took about two months to write and, apart from design- 


ee 
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ing an optimum system, it has reduced the actual design 
time from 3 days to about 25 minutes. 

Input to the program consists of two types: 

© Constants. These are the cost values and allowable 
minimum and maximum velocities. They can be updated 
from time to time, but not necessarily for each program. 

© Regular input data for the program. These consist 
of description of lines through the use of code numbers 
already assigned for various fittings, the description of 
circuits and their lines. 

The input and output shall not be described in detail; 
however, the printout of an example problem is pre- 
sented. The whole thing is self-explanatory except for 
input data that uses already assigned code numbers for 
the description of the fittings in the line. 


Example Problem. A cooling water system for a plant 
at Mims, Fla., was designed by the computer program. 
The cooling water system consists of 42 lines and 11 cir- 
cuits. The available pump head is 108 feet and static 
head is 20 feet. The program took about 25 minutes to 
design the system. The input data and the results of the 
program are presented in the accompanying tables. 

The program does the optimum sizing of the system 
based on installed cost. The results of the example prob- 
lem show that the various circuits are very nearly bal- 
anced and are very close to the allowable limit, Eight out 
of 11 circuits are within 10 percent of the allowable 
value. Also, the lines are sized such that the system cost 
is the least. This is the most significant aspect of the 
program, for in sizing the line the program takes into 
account every fitting present in the line and optimizes 
the over-all cost. 


Advantages of the Program. 

© Initial selection of the pump for a given system. A 
routing of piping can be run for various allowable heads 
and the cost of the systems can help in selection of the 
pump. 

© Design of the system 

© Calculation of the cost of the systems. 

The program can be further extended to prepare the 
field and fabrication bill of materials, purchase requisi- 
tions, and isometric drawings. Using the same input data, 
all of the above items can be obtained. Work is presently 
being done in the above field and very soon the pro- 
gram shall be able to prepare the field and fabrication 
bill of materials. 

The program makes full use of the advantages offered 
by the computer over the hand calculations. It saves in 
the computational time as well as in the cost of system. 
It completely eliminates the guess work that is usually 
involved in the sizing of a system by hand calculations. 
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Simplified Utility Loop Balancing 


Utility distribution systems require 
periodic checks on flow and pressure 
drops in the loop mains. This 
method reduces the calculation time 
from days to minutes 


B. West and A. J. Newton 
Polymer Corp., Ltd., Sarnia, Canada 


Urmiry pIsTRIBUTION sysTEMs for refineries and chem- 
ical plants are becoming larger and more complex. Usually 
these systems are looped, so that the utility supply to any 
unit is not dependent on one path through the system. 
As a result, complicated piping networks are built up. 

A major plant expansion will normally include a study 
of the utility distribution systems, to determine the effect 
of the load increase, and to indicate where expansion is 
required, These studies when performed by hand are 
very time consuming and prone to error because of the 
large number of repetitive steps involved. 

The repetitive nature of the calculations naturally leads 
to the use of a computer method of evaluation. Com- 
puter programs 2, 3, 6, 7, 8 have been written to solve 
this type of problem. 

The following is a description of the network analysis 
program developed at Polymer Corp. 


Information Required. It is decided to examine one of 
the plant utility distribution systems. The object of the 
study is to determine the effectiveness of the system under 
present and future loads. The effectiveness might be de- 
termined by comparing required pressures at various 
points in the system to the predicted pressures. It could 
also be determined by comparing the predicted fluid 
velocity in the pipes to the maximum allowable fluid 
velocity. 

This information can be provided by determining the 
flows and pressure drops in the pipes of the system. The 
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procedure used to find the correct flows is relatively 
simple but tedious to solve by hand. The calculation pro- 
cedure can be broken up into steps. 


1. Collect Data. Pipe diameters, pipe equivalent lengths 
(ie. including length for fittings.), system configuration 
and average fluid properties. 


2. Estimate Flows in the Network. The flow rate of 
fluid in all the pipes in the system is estimated. The only 
criterion for this estimate is that the chosen flows must 
be consistent with the over-all mass balance. It is pref- 
erable, however, to make reasonable estimates, as this will 
reduce the length of the calculation. 


3. Calculate the Pressure Drop in the Pipes. Any 
standard pressure drop calculation is used. In the com- 
puter program, the Polymer standard procedure is used. 
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END 


Fig. 1—Simplified computer flow diagram. 
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4. Check to Find if the Loops Balance. Each loop is 
examined individually. The algebraic sum of the pressure 
drops in the loop pipes is found. If this is zero, the loop 
is in balance. If this sum is not zero, an adjustment must 
be made to correct the flow rates in the loop, to give 
balance. Each loop is checked and corrections made to 
the estimated flows if necessary. 


5. Return to Step 3. Steps 3 and 4 are repeated until 
all the loops in the system balance. 


The above procedure is applicable to both hand or 
computer solution. Hand calculations make use of various 
approximations because a detailed calculation would take 
too long. Computer solutions may be made rigorous by 
reducing the number of approximations. 


BASIS FOR COMPUTER CALCULATIONS 


The relationships used for pressure drop calculations 
in our program are simplified and reflect certain assump- 
tions stated below. The degree of sophistication used to 
calculate the pressure drops is very much dependent 
upon the individuals’ needs. For our requirements we 
have found the flow formulas presented here to be ade- 
quate. 

The assumptions inherent in the following theory are: 


@ Isothermal flow for compressible fluids, 


© Fluid properties are for the average temperature and 
pressure of the fluid in the system, 


© The pipes are rough, and 
© Datum temperature is 50° F. 
The basic relationships used for analyzing the loops are 
the equations for flow of fiuids in pipes. 
Pressure Drop Equations 


Noncompressible Fluids. The relationship between inlet 
and outlet pressures is 


Le \ ( G? 
ne 2)(2) us ae 


OF fy — by = ky (S.G.) fLy (=) (2) 


) 








— 9 =4.7 (c8)(V0) (£2) Q (4) 
Equation 2 and 4 can be put in the form: 


H,=Ki; cay ($8 ) (5) 


The value of K will be constant for a given system. 
Equation 5 is the general equation used in the com- 
puter program to calculate pressure drops in the pipes. 


Friction Factor. In Equation 5 all the terms except f, 
the friction factor, are input data predetermined for the 
system under study. The friction factors can be read off 
a chart if the calculation is done by hand. The computer 
solution requires an analytical expression for friction 
factor. An expression developed by Dukler* from Moody’s 
friction factor chart and used in the program is shown: 


ap = 108 + 5.76 og ae 
ac 2 2 
— 1.75 log R — 1.10 (log R)? (6) 
=Re () a 
D 2 
Substituting the following into Equation 6. 
f=4f 
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where 


gives: 
= = 1.03 + 5.76 log W — 1.75 log R — 1.10 (log R)? 


or 
FY = 1.03 ++ 2.5 In W — 0.76 In R — 0.207 (In R)? 
(7) 


The correct estimate of f gives FY = O. 

To solve Equation 7, a value of Y is chosen and FY 
calculated. If FY is not zero, a correction is made to the 
initial value of Y. This correction is found by application 
of the Newton-Ralphson technique.’ 








Compressible Fluids. 
Where 
G Le 
bf —h2= ate EE (4 :)] (3) ¥,=¥,.,-E 
2 Py and 
1 
_ FY _ 103—\¥ ) 42.50 (In W) — 0.76 (In R) — 0.207 (In R)? 
~~ d(ER =e 1 1.74 0.414 
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For long pipes (-+ 200 ft) or low 


velocities In (=) approaches zero. 
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This procedure is repeated until FY is zero then f is 
found from Y. 


Balancing Loop Flows. The pressure drops in all the 
Pipes are calculated based on the initial estimates of flow 
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rates, The algebraic sum of the pressure drops is found 
for each closed path or loop taking the clockwise flows as 
positive (this is the Hardy-Cross convention) . If the esti- 
mated flow rates are correct, this sum will be zero for all 
the loops. If this sum is not zero, an adjustment must be 
made to the estimated flows. The adjustment is made ac- 
cording to the Hardy-Cross® formula which is a special 
application of the Newton-Ralphson technique. 


isin icine . 
Algebraic Sum = ASj= )) Hy = ) [ Atwfen te | 
ist i=ii ii 


(9) 


If AS is not equal to zero a correction is made to the 
estimated flows in the loop, j. 
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The new estimated flows are used to recalculate the 
pressure drops and the procedure continues until the loops 
are balanced. 


ie. until 


j=m 
y | 4s;|=0 
ja 


Tolerances on the Iterative Procedures. Iterative 
techniques are used to determine the friction factors and 
the correct flows in the network. The criterion for finding 
the correct value is that a test variable (FY and & AS;) 
should become zero. An optimum tolerance is required 
on the convergence test of these variables. A balance must 
be struck between the accuracy of the results and the 
computing time. 


Friction Factor Calculation. Two different tolerances 
are used on the test variable FY. Initially a tolerance of 
xt 0.1 is used and if the flow rates are changed by less 
than 5 percent from the flows in the previous iteration, a 
new friction factor is not calculated. When the flows are 
almost balanced the friction factor criterion is tightened. 
The tolerance becomes =t 0.001 and a new friction factor 
is calculated for each change of flow rate. 


Flow Balancing Criterion. Ideally, when the loops are 
in balance the values of AS in each loop equals zero. To 
attain this value would require considerable computing 
time. The criterion used to decide whether the flow rates 
need further modification is to sum the absolute values 
of the deviations from zero of each loop, and if this 
exceeds a desired value, a further change is made to the 
flow rates. This desired value is best determined anew 
for each system, taking into account the computing time 
available and the desired accuracy of the results. 


Absolute Pressure Calculations. The pressure drops 
in each main in the network are now known. For com- 
pleteness, it is necessary to calculate the absolute pressure 
at the end of each pipe. 
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The following systematic approach is used. The loops 
are treated in turn, Loop one contains the input (or the 
highest pressure input if there are several inputs) to the 
system, at which point the pressure is known. All the 
pressures are calculated around loop one. The next loop 
in numerical order is then considered. A pipe is found 
which is common to a previous loop or enters an intersec- 
tion common to a previous loop. This search provides a 
point of known pressure from which the other pressures 
around the loop may be determined. All the loops are 
examined to give all the pressures throughout the system. 
It is important to number the loops in such a way that 
they have mains or intersections in common with lower 
numbered loops. This is necessary to ensure that each 
loop has a starting point where the pressure has already 
been found, 

Fig. 1 shows a simplified flow diagram of the computer 
program. 


Experience With the Program. The program described 
here has been used extensively by the Utilities Department 
of Polymer Corp. for evaluating proposed changes to the 
steam and water distribution systems. The program has 
been checked by comparing the results obtained with the 
actual system. A good agreement has been obtained. 

The great reduction in calculating time has enabled 
more varied modifications to be studied for each expansion. 

A sample problem follows which shows the preparation 
and computing times for the problem and for typical 
utility problems. For the example shown, the computer 
gives the result in about 15 minutes; whereas the hand 








SIMPLIFIED UTILITY LOOP BALANCING . . . 


2,000 U.S. gpm 





1,000 U.S. gpm 3,000 U.S. gpm 


10,000°U.S. gpm 4,000 U.S. gpm 


Fig. 2—Sample problem. 


TABLE 1—Input Data for the Sample Problem 











Main Number Diam. (in.) Eq. Lath. (ft.) Est. Flow USgpm 
1 12 1000 5000 
2 10 600 —5000 
3 10 500 4000 
+ 10 1200 —5000 
5 10 1500 —3000 
6 10 200 4000 
7 10 1000 1000 





TABLE 2—Computer Calculated Values of Flow and 
Pressure in Each Loop Main 








Main Number Balanced Flow, USgpm —_ Outlet Pressure, psia 
1 5934 68.2 
2 —4065 80.7 
3 4934 34.6 
4 —4065 34.6 
5 —2235 16.4 
6 4764 24.1 
iz 1764 16.4 





calculation would take at least one hour, The advantages 
of the computer are greatly increased as the size of the 
system is increased. One hour’s computing time replaces 
two weeks or more required for hand calculations on the 
big utility systems. 


Sample Problem. Determine the flow rates and pressures 
in the network shown on Fig. 2. 

Table 1 gives the main numbers in the loop and tabu- 
lates the diameter of the main, the equivalent length, and 
the estimated flow in USgpm. 

Table 2 is the computer calculated flow in each main 
and the outlet pressure. 

The time taken to complete the problem was as follows: 


Time 


Function 


Prepare data 
Punch data 
Computer solution 


10 minutes 
2 minutes 
3 minutes 


15 minutes 


Hand vs Computer Calculation Time. A steam dis- 
tribution system consists of 29 mains and 7 loops. 

A typical computing time for this system based on a 
reasonably accurate initial estimate is 30 minutes. 
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The service water distribution system consists of 41 
mains and 10 loops. 

A typical computing time is 15 minutes. 

No specific times of calculation by hand are available, 
but several days of continuous work were required to pro- 
duce results of comparable accuracy. 

A recent paper® gives an excellent description of a much 
more rigorous method of solution to complex pipe system 
problems. The method has obvious attractions for use 
in systems where the fluid conditions and properties are 
subject to considerable change. The simple procedure 
presented here, appears sufficiently accurate for fluids at 
steady conditions. 
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NOMENCLATURE 


AS Algebraic sum of pressure drops in a loop 
CF Compressibility factor of fluid 
D Diameter of pipe, inches 
f Weisbach-Darcy friction factor 
f} Fanning friction factor 
FY Friction factor test variable 
g Newtons-law conversion factor, ft.-Ibs. 
force-sec.* 
G Mass velocity, Ibs. mass/ft. sec. 
H_ Pressure drop in pipe, lbs. force/in.? (Ibs.-force?/ 
in.4) (units for compressible fluids) 
k, Constant for noncompressible fluids, Ibs. force in.> 
min.?/ft. USG? 
k, Constant for compressible fluids, 
hr.*/ft.8 °R Ib, mass 
K Constant for system and fluid 
K1 Absolute roughness magnitude 
Le Equivalent length of pipe, ft. 
M Total number of loops 
N Total number of pipes 
, Inlet pressure, Ibs. force/ft.2 
. Outlet pressure, lbs. force/ft.* 
Q Volume or mass flow rate, USgpm or (Ibs. 
mass/hr.) 
Re Reynolds no. 
S.G. Specific gravity of fluid at flow temperature 
T Absolute temperature, (°R) 
VO Specific volume of gas at standard temp. and press. 
(ft.3/Ib. mass) 
p Fluid density at flow temperature, Ibs. mass/ft. 
4Q_ Flow correction factor, USgpm or (Ibs. mass/hr.) 
E Friction factor correction factor 


mass/Ib. 


Ib. force? in.® 





SUBSCRIPTS 


i. Pipe number 

j Loop number 

m Over-all iteration number 

I Friction factor iteration number 
jn Number of pipes in loop j. 
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Figure Optimum Absorption Design 


The minimum cost to recover solute 
can be determined graphically from the 
computer data plotted here 


4. C. Kim and M. C. Molstad 
University of Pennsylvania, Philadelphia 


Tue minimum unrr cost of recovered solute is ob- 
tained by the proper selection of absorber height and op- 
erating conditions. The relationship of these optimum 
conditions is presented graphically. To use the method, 
it is necessary to evaluate the quantity C; © r m Gu/Cs 
Hoo (Kn —1). An example is solved showing the use of 
water to recover acetone vapor at minimum unit cost 
from an air stream. 


Governing Equations. Equations giving the total annual 
cost of treating a gas stream entering an absorber have 
long been available’ and have more recently been put 
into somewhat modified form.** The total cost includes: 


© the fixed cost of the absorber, to which has been 
added for convenience one-half of this fixed cost to cover 
the energy required to overcome gas pressure drop 


© the value of the unrecovered solute taken at the 
market cost of its replacement 


© the cost of recovery of the solute in the stripper op- 
eration. 


The cost to be minimized is the total annual cost of 
recovering one pound mole per hour of solute vapor fed 
to the absorber in the gas stream, i.e. Cr/GuS (Y:— Y:). 
The solute content is assumed to be small enough so that 
the total molal gas flow rate, Gy, can be considered con- 
stant in the absorption column. The expression for the 
foregoing cost is 





Cy ee Cr 
GuS i= 92) ~~ GySy, 1—1/¥) 
CZ ‘Gi 





= Gpattaiyy) =) 


r{1—[(KpGy/Ly) (91) (1—1/¥)1} +1) 


C, a eee 
+ C:9\ Tee 1) (Gu/Eu) 01) G17) 


(1) 
where Z is the height of absorption tower given by 


In { [1 — (mGy/Ly)1 ¥ + (mGy/Lu)} 
1— (mGy/Ly) 


Equation (1) is partially differentiated with respect to 


Z=Hog (2) 
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Fig. 1—A digital computer was used to establish these curves. 


Wak th 










$0 t 
Yen ath 
Fig. 2—Results are plotted to show the minimum design. 


the absorption factor, and the resultant equation is set 
equal to zero. The final expression, from which the opti- 
mum absorption factor is computed for the foregoing case, 
then becomes: 
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(Ly/mGy) ((Ly/mGy) — 1] = 


ie 
[C,Hoq (Kp — 1) /C, @rmGy] ( 


A Direct Solution. All the quantities in the coefficient, 
the first term on the right side of Equation (3), are given 
or known in an absorption operation, The remaining 
terms are explicit functions of the absorption factor 
Ly/mGy and Y, the ratio of the inlet and exit solute con- 
centrations in the gas stream. Therefore this equation can 
be solved for an optimum absorption factor as a function 
of Y and the coefficient. Previously this has required a 
series of cut-and-try calculations. This paper presents a 
direct solution. 

The actual computation was carried out by an IBM 
7040 digital computer, and the results are presented in 
Fig. 1. For convenience the figure shows the inverse of 
absorption factor versus the inverse of the coefficient, with 
Y as a parameter. The related pairs of values of the ab- 
sorption factor and Y, all satisfying Equation (3) for any 
lean absorption, may now be used in Equation (1) to 
locate the minimum. 


An Example. The illustration of the removal of acetone 
vapor from an air stream by water in an absorber packed 
with 1-inch Raschig rings’? can be used. The following 
same values were used: for the acetone-air-water system 
m = 2.7 and Ky = 23; a mass velocity of 735 pounds of 
gas/(hr.) (sq. ft.), equal to 25.4 pounds moles of air/(hr.) 
(sq. ft.) and a corresponding Hog of 2.5 feet; r = 1,25 
in the stripper; © = 8,400 hours per year; and C; 
$4.35 per cubic foot of absorber per year. Values of 
$0,018 per pound-mole of steam produced at the base 
of the stripper and $4.90 per pound-mole of acetone were 
used as today’s costs. The value of 0.02 taken for y, was 
chosen to be less than the flammability limit of acetone- 
air mixtures, 

The inverse of the coefficient in Equation (3) equals 
54.1, and Fig. 1 gives at this abscissa the corresponding 
pairs of values of Ly/mGy and ¥ which satisfy Equation 
(3). A few of these pairs are substituted in Equation (1) 
and the results plotted in Fig. 2 to locate the minimum, 
This is at ¥ = 200, with a total annual cost of $1,865 for 
the recovery of acetone at the rate of one pound-mole per 
hour of operation. 





Interesting observations can be made from these results, 
The cost for the recovery of acetone is 


1,865 ($/yr.) /lb.-mole acetone) /(hr.) 
8,400 (hr./yr.) x 58 (Ib.)/(l 


= $0.0038/Ib. acetone 








-mole acetone) 





or less than one-twentieth of the current market price 
used in the example. 

The values for the three cost terms in Equation (1) at 
optimum conditions are $430, $205 and $1,230, respec- 
tively. At Ly/mGy = 1.22 the height of the absorber is 
20 transfer units or 50 feet. When ¥Y = 200, 99.5 percent 
recovery of the entering acetone is obtained. (The final 
recommended design would undoubtedly be for a shorter 
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[1 — (mGy/Ly)] ¥ + (mGy/Ly) 


In {1 — (mGy/Ly)1 ¥ + pie) (3) 


1— (mGy/Ly) 
tower, where the incremental return on the added invest- 
ment would be the minimum acceptable.) 

The assumed values of m and Gy give a water mass 
velocity of 1,510 pounds water/(hr.) (sq.ft.), which is 
found to, at or near a loading condition,‘ although the 
pressure drop is only 0.8 inches of water per foot of tower 
packing. Therefore the calculation should be repeated 
with a lower Gy. 

For any problem the first cost term is affected by the 
values taken for Hog and Cs. In the absence of experi- 
mental data the published Hog values for the dilute 
ammonia-water system could serve as a guide.> The value 
of C; would be taken from current costs of absorption 
equipment.®715:8 


Thus the optimization of an absorption operation is 
determined by evaluating the various terms in the dimen- 
sionless quantity C,@rmGu/CsHoo (Ky — 1). For any 
value of this quantity the height of the absorber and its 
operating conditions for minimum unit cost of a recov- 
ered solute can be readily determined. 


Nomenclature 
C, Annual cost of power and fixed charges per volume 
of absorption column, $/(cu. ft.) (year) 
C, Cost of solute, taken at current market, $/Ib.-mole. 


C, Total cost of stripping operation, expressed as 
$/lb.-mole of vapor produced at the base at the 
stripper; includes fixed charges, cost of cooling 
water, and cost of steam 

C, Total annual costs, $/yr. 

Gy Molal gas velocity through absorber, 1b.-mole/ 
(hr.) (sq. ft.) 

Hog Height per transfer unit, ft.; assumed constant 
Slope of the equilibrium curve (y= Kx) in the 
stripping column; assumed constant 
Molal absorbent velocity through absorber, 1b.- 
mole/(hr.) (sq. ft.) 
m Slope of the equilibrium curve (y= mx) in the 
absorption column; assumed constant 
mGy/Ly Absorption factor 
r Ratio of actual to minimum reflux ratios in the 
stripping column 
s Absorber cross-sectional area, sq. ft. 
* Mole fraction of solute in the solvent 
Yys¥ Mole fraction of solute in the gas stream at the 
inlet and exit of the absorber, respectively 
I/ Vo 
Hours of operation per year 


on 
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Best Approach to Compressor Performance 


When a centrifugal compressor is 
specified for nonideal gases 

and the manufacturer’s data is based 
on air, how do you predict design 
and off-design performance? 

This calculation procedure solves 

the problem 


Stephen A. Shain 
Shell Development Co., Emeryville, Calif. 


THE sPEcIFICATION of a centrifugal compressor and 
prediction of design and off-design performance are parts 
of a problem that continually arise in process design 
studies. To make these calculations less difficult, a com- 
puter program was written which handles very general 
situations including multistage compression, inter- and 
intra-stage leakage, multicomponent gas streams, depar- 
tures from perfect gas behavior, intercooling when nec- 
essary, and allowances for manufacturers’ data on the 
basis of air to permit design for hydrocarbons and other 
vapors. 

The equations used for the individual compression 
stages of the compressor are summarized in the sections 
which follow. The assumptions and limitations implied 
in the general efficiency and head coefficient correlations 
used, are discussed. The equations given are more general 
than those usually appearing in the literature.t)*%®%1* 


COMPRESSOR PERFORMANCE EQUATIONS 


Outlet Temperature and Pressure. The temperature 
and pressure at the outlet of the compression stage are 
related to the hydraulic efficiency (polytropic efficiency) 
and head coefficient by the following set of equations: 


(0.0012854) A @inZ 
a (60)?G, Le + \omT), 
7,\" 
P,=P,(>2 (2) 
1 


where 





T, =T, 








100 


The efficiency and head coefficient* are defined by 
« — (0.0012854) Wy 


a= (4) 
2 (60)? 2, W, 
eae eee 
72D2N2 (5) 
where the polytropic head, Wp, is defined 
We= aw [2 (Polytropic path) (6) 
p 


The properties appearing in Equations (1) and (3) 
are to be averaged over the compression path. The tem- 
perature equation, Equation (1), requires an efficiency 
which is an appropriate average over the path of the 
compression, but the pressure equation is only valid for 
a constant efficiency path. 


Dimensionsless Numbers. Dimensional analysis of the 
compressible fluid flow equations’ indicate those variables 
and groups of variables which are important in specify- 
ing the performance of a centrifugal compression stage. 
Applicable equations are the equation of of continuity, 
the momentum equation, the energy balance, and the 
equation of state. These equations must be solved subject 
to specified boundary conditions. 

The solution to the equations is controlled by the 
magnitude of the following dimensionless groups: 


The Grashof number 


gBATL! 


er 


The Reynolds number 











(8) 

The Prandtl number 
(9) 

The Eckert number 
a (0.012854) U2 (10) 


aC, AT 


If the reference temperature, AT, in the energy bal- 
ance equation is taken to be the absolute temperature, 
T, the Eckert number is related to the Mach number 


b+ Ge] 


_ * The definition of the head coefficient, Equation (5), involves a factor of 
2 which is omitted in some sources. 
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Fig. 1—Head correlations for a variety of speeds. 


In the above list of dimensionless groups, the Eckert 
number (Equation 10), may be replaced by 


The Mach number 
U 

Nua = We (12) 
The speed of sound, appearing in the denominator of 
Equation (12), is the speed at which a small amplitude 
disturbance will propogate through the gas in the ab- 
sence of absorption losses, i.e., for an adiabatic reversible 
or isentropic process. The speed of sound is given by the 


relation: 
oP 
=| 144 
vite (7), 


The relations of thermodynamics may be used to 
transform this equation to a form which allows direct 
use of the equation of state: 


i— =~ 
V,= ite) 
8e}r 


The final relationship for the speed of sound in a 


nonideal gas 
v= 144 ge KZRT 
a 
esl — InZ 
ainP }, 


The calculation makes use of Equation (12c) evalu- 
ated at the inlet of a compression stage. 


v 





(12a) 








(120) 





(12c) 





In the computer program, the characteristic velocity 
(U) is taken to be the impeller tip speed. In addition 
to the above groups, the equation of state and the 


boundary conditions are implicit in determining the 
solution, 
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Two compression stages will have the same solution 
in terms of appropriate dimensionless variables provided 
certain side conditions are satisfied. First of all, the 
boundary conditions in the two situations must be iden- 
tical. Second, both compression stages must be geometri- 
cally similar, that is, the ratios of all dimensions must 
be identical, and third, the equation of state for the 
two gases undergoing compression must be identical. 


Thermodynamic Variables. Solutions to gas motion 
equations are generally given in values of gas thermo- 
dynamic properties as functions of spatial variables 
in the impeller and diffuser. However, the solutions of 
interest are certain functions of these thermodynamic 
variables rather than their complete spatial variation. 
In particular, the solution may be thought of as the 
ratio of useful work to the total work or total enthalpy 
rise, and the ratio of useful work or pressure head to the 
kinetic head of the impeller. This first variable is the 
polytropic head divided by the actual work input and is 
referred to as the hydraulic efficiency. The second vari- 
able is the head coefficient. These variables have been 
previously defined by Equations (4) and (5). These ratios 
will depend on the boundary conditions and in particu- 
lar on the ratio of the fluid velocity to the tip velocity of 
the impeller. This ratio is the volumetric flow coeffi- 


cient (4) 


Q 


a@DNA (3) 


and the formal relations for the efficiency and head 
coefficient and the flow coefficients are: 


e= 8 ($;NopNpeNpp Nya equation of state 
parameters, geometric variables) 


Y= ($;Nop Nee New Nyai equation of state 
parameters, geometric variables) 


The volumetric flow used in defining the flow coef- 
ficient may be an inlet flow, an outlet flow, or any aver- 
age of inlet or outlet flows. Several alternate approaches 
have been reported in the literature.***'%"! Better cor- 
relations appear to result from a use of intermediate 
flows, but correlations in this form are not always avail- 
able from the manufacturers. The example given is 
based on inlet volumetric flows, 

For impellers which are geometrically similar, the 
area (A) is uniquely related to the impeller diameter 
(D). However, the value of A has been left free so that 
it may be adjusted to improve the efficiency and head 
coefficient correlations or to put them in a convenient 
form, particularly if geometric similitude is not main- 
tained. For example, the flow coefficient may be inter- 
preted as the ratio of the volumetric flow to the manu- 
facturers’ design volumetric flow if a nominal area (A) is 
taken to be 
Q (design point) 


A= 
aDN 


(14) 
This form is often convenient since the manufacturers’ 
correlations are commonly in terms of percentage flows 
relative to the design point, i.e., Equation (14) fixes ¢ = 1 
at the design point. A similar definition of a nominal 
diameter (D) may be obtained from Equation (5) for 
wv =1 at the design point. 
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Efficiency and Head Coefficient Correlation. When 
the heat loss to the surroundings is negligibly small, the 
Prandtl number may be dropped from the solution.” If the 
Reynolds number is large (Nee* >> Ner), natural con- 
vection is no longer of importance; thus the Grashof num- 
ber may be eliminated from the solution.” If the Mach 
number is smaller than 1, that is, if the speed of sound 
is not exceeded, at all points along the compression path 
the performance of the compression stage is insensitive 
to the Mach number.®"* For Reynolds numbers greater 
than 100,000, as is usual, the results are insensitive to 
the Reynolds number.**® 

The mechanical design of different impellers for a 
given frame or frames is not generally based on the 
principle of geometric similitude. Impellers for the same 
frame designed for varying volumetric flow capacities 
will typically have the same impeller diameter with a 
varying channel width or area of flow. The assumption 
of geometric similitude is then not strictly applicable. 
‘An assumption implied in the example calculation is 
that a single correlation for head coefficient and for 
efficiency applies for all of the impellers available in 
the computation. This assumption will be strictly valid 
only if the impellers are geometrically similar. However, 
it is necessary to make this assumption in many cases 
since more detailed information on the performance of 
the individual impellers is not always available. 

The form of the correlations which is suggested for 
use reduces to: 


e=e (9) 
¥=¥ (9) (15) 


An example correlation of head coefficient as a func- 
tion of volume coefficient is presented based on the data 
of C. A. Macaluso‘ for a single stage compressor. Since 
the actual physical dimensions of the compressor were 
not included as part of the data, the results are corre- 
lated only in terms of ratios of the appropriate quanti- 
ties. That is, the head coefficient is replaced by the 
dimensional ratio (in enthalpy units) 


H Btu/lb mass 
raw ke Pevbmnae, 
nN [ (rpm)? ] cs) 


and the volume coefficient (based on suction volume) is 
replaced by the dimensional ratio 
_ Q, emf 


fe 'N 7pm 





(17) 


The results for a variety of speeds ranging from 5,400 
rpm to 8,250 rpm and volumetric flows of 2,000 cfm to 
4,500 cfm are presented in Figure 1. While there is scat- 
ter of the results, a single correlation can be presented 
with an accuracy of 5-10 percent. A point of particular 
interest is that the results are for compression of differ- 
ent gases (molecular weights of 120.9, 137.4, and 170.9). 
The surge limit or minimum value of the volumetric 
flow coefficient and the choking limit or maximum value 
of the volumetric flow coefficient are reduced to single 
values of the volumetric flow coefficient when presented 
on the basis of Figure 1. 

The above results and others from the literature*** 
suggest that correlations for both head coefficient and 
efficiency as functions of volumetric flow coefficient may 
be used for a variety of nonideal gases provided the ap- 
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Fig. 2—Block diagram of compressor stage 























Begin : 
a calculations. 
Initialize 
Outlet 
Stream 
Calculate <max 7 Count _>rnae 
Outlet Iterations >— 
Properties 
ae ox >0 | Calculate heat 
est 0 Teste Saree losure 
<0 <0 
rece Store for 
Next Stage 














propriate dimensionless forms are retained. In the ab- 
sence of information on the particular gas being com- 
pressed, manufacturers’ information based on air can 
be used in the computation even though exact similitude 
would not be obtained. 


Compressor Performance Calculations. The perform- 
ance calculation for a compressor or series of compressors 
is accomplished on a wheel-by-wheel basis. The perform- 
ance of each impeller or compression stage is considered 
in sequence by application of the general set of perform- 
ance equations, Equations (1), (2), and (3). These equa- 
tions allow calculation of the temperature and pressure 
at the outlet of a compression stage provided that the 
efficiency and head coefficient are known. 

The specific definitions of the coefficients and dimen- 
sionless groups used have been presented as Equations 
(4) and (5). The efficiency and head coefficient are 
represented as functions of a single variable, the volu- 
metric flow coefficient, Equation (15). The correlation 
of efficiency and head coefficient with volumetric flow 
coefficient is often a satisfactory representation of the 
performance of not only a specific impeller, but also of 
other impellers of similar design but which are scaled 
for use at different volumetric flows. This correlation of 
efficiency and head coefficient provides a convenient 
condensation of information on the performance of a 
compression stage for a variety of operating conditions. 

There are limitations inherent in this approach, and 
some of these have been discussed above. The calculation 
procedure is sufficiently flexible so that additional in- 
formation, when available, can be included. 


Outline of Program Logic. The program logic may be 
schematically outlined as follows: 


1. Read in input data. 
2. Calculate recycle (leakage) flows. 


[If the compressor calculations are incorporated with 
other process calculations, recycle flows (not necessarily 
leakage, as in this example) may depend on the com- 
pressor discharge conditions. In that case an over-all 
material balance may require iterative application of the 
entire compressor calculation. ] 











10 2 ® \ 
Stage ‘Stage Stage 
ees re 20s: 7 t 




















Fig. 3—Flow diagram for a three-stage compressor. 


3. Calculate physical properties for the suction. 
4. Initialize stage discharge physical properties. 

5. Calculate appropriate average flow coefficient. 

6. Calculate head coefficient and efficiency. 

7. Calculate new discharge temperature and pressure. 


8.Compare current and previous discharge temper- 
ature. If closure to within a specified tolerance is 
not achieved, return to (5). If closure is achieved, 
calculate the next stage starting with (4). After 
the last stage, continue with (9). 


9, Print out final results. 


Figure 2 is a simplified diagram of the subroutine 
which calculates the performance of a single stage. 


Use of the Computer Program. The input data re- 
quired by the program can be divided into the following 
categories: 


Component data—consisting of a specification for 
vapor components in the form appropriate to the sub- 
routines which calculate the compressibility factor, Z, 
and other physical properties required for the compu- 
tations. 

The equation of state used in the example is based 
on a modified version of the Martin and Hou equation 
of state’ 1° 





Ay + BT + CyeK?/7, A 
(= 5)s + GH 





Impeller correlation data—specifying a correlation of 
efficiency and head coefficient of the impeller or impel- 
lers with the volumetric flow coefficient. (The correla- 
tion may be based on inlet flow, outlet flow, or any 
average of inlet and outlet flows from a compression 
stage. Variations in the efficiency and head coefficient 
with Mach number and Reynolds number may be in- 
cluded, if desired.) 


Compressor sizing data—consisting of a list of impeller 
scale parameters or physical dimensions to be used in 
the correlations for efficiency and head coefficient. 
(These usually are a characteristic impeller area and 
diameter, but other variables which reflect the capacity 
of the impeller may be used.) 


Stream data—consisting of a specification of the inlet 
stream temperature, pressure, and molar flows for each 
component. 
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Miscellaneous specifications—including the rotational 
speed for each frame, and several bounds and increments 
which must be specified. 


The output from the program is a reiteration of the 
input data and a statement of the calculated properties 
for the inlet stream, i.e., the density and the volumetric 
flow. A listing of parameters and properties for the out- 
let stream of each compression stage is presented. These 
include the stream number, the frame number and im- 
peller number, the nominal diameter and area for the im- 
peller, the rotational speed, the temperature, pressure, 
density, molecular weight, the volumetric flow, and 
arithmetic average value of the specific heat ratio, the 
value of the exponent which defines the compression 
path (given by Equation 3), the value of the average 
flow coefficient, the Mach number, Reynolds number, 
head coefficient, efficiency, the head and horsepower. 
When a stage consists of intercooling, the temperature 
and pressure of the stream is presented, and, if phase 
separation occurs, molar flows of both the vapor and 
the liquid stream are presented. 


In addition to the calculated properties, a warning 
flag, FLAG, is indicated for each compression stage. 
Values of FLAG are used to indicate the course of the 
computation and signal any difficulties which may have 
been encountered during the computation. The warn- 
ings which may be indicated by FLAG are: 


© The Mach number exceeds the critical Mach number, 


© The average value of the flow coefficient is less than a 
minimum value. 


© The average value of the flow coefficient is larger than 
a maximum value. 


© The number of iterations for the calculation at this 
stage has exceeded the maximum number of iterations, 
i.e., the calculation has failed to converge. 


© The calculated head coefficient is negative (may result 
from inappropriate extrapolation of correlation). 


© The calculated efficiency is negative (as in head coef- 
ficient above). 


© Intercooling has resulted in condensation. 


© The temperature has exceeded the maximum allow- 
able temperature. 


Example Calculation. An example is presented on the 
application of the above procedure using a large scale 
digital computer. Calculations are made for a three-stage 
propylene refrigeration compressor. The shape of the 
head-capacity and the efficiency-capacity relations was 
taken to be the same for all three wheels, with the flow 
coefficient based on inlet flows. The head-capacity rela- 
tion for the second stage was used as supplied by the 
manufacturer. The shape of the curves for the first and 
third stages differ by 2-5 percent at the extremes in the 
range of 60-140 percent design capacity. An efficiency- 
capacity relation was not available for the wheels in 
this machine. The manufacturer was able to supply an 
efficiency-capacity curve for a similar machine at a 
specific speed which was close to the design condition 
at the second stage. This curve was used to calculate 
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Table 1—Example Calculations Compared With Performance Test 
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‘Temp., 
°F psia cfm> 


Flows 
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DISCHARGE (3rd Stage) PERFORMANCE TEST 
Press., 


Brake, Brake, 
hpe psia hp 


psia 








#-100 Percent Propylene. 


an efficiency-capacity relation for the range in capacity 
of from 60-140 percent of the design points. Again, the 
wheels may differ by approximately 5 percent at the 
extremes of the range. The flows in the recycle (leak- 
age) streams were taken as 1.5 percent at each stage and 
6.5 percent from the discharge of stage 3 to the suction of 
stage 1. Figure 3 is a schematic diagram indicating the 
interstage flow, stream 12, and the intrastage flows, 
streams 8, 9, and 10. 

The results of the computations are summarized in 
Table 1 and are compared with a performance test. 

The discharge pressures of 312-276 psia and power 
requirements of 2949-3107 hp are somewhat above the 
275-200 psia and 2,400-2800 hp ranges from the original 
design specifications. However, the calculated values are 
close to those from the actual performance tests for this 
machine. The calculated discharge pressures (third 
stage) are within 1-3 percent of the pressures actually 
achieved with even closer agreement for the brake horse- 
power. The agreement is excellent considering that man- 
ufacturers’ head and efficiency data for air has been 
used to calculate performance for propylene under con- 
ditions where it is a nonideal gas. 


Why This Approach? The approach presented here is 
intended as a thermodynamically consistent method for 
performance calculations when nonideal gas phases are 
important. It is hoped that the presentation will en- 
courage publication of further experimental results which 
will allow evaluation of this approach and outline more 
clearly its limitations. The author would like to empha- 
size to those concerned with design, manufacture, and 
performance of such compressors that methods and com- 
puter capabilities for evaluation of thermodynamic prop- 
erties of fluids are becoming generally available and 
should extend the usefulness of many design methods. 
The physical properties programs'® used here have been 
made available to the Physical Properties Project Sub- 
committee of the American Institute of Chemical Engi- 
neers. 
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>For feed flow, temperature, and pressure. 


=1.02 x Compression hp. 
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NOMENCLATURE 


nominal area, sq. ft. 

heat capacity at constant pressure, Btu/Ib, mass-°R 
nominal diameter, ft. 

local acceleration of gravity, ft./sec. sq. 
gravitational constant, 32.16 Ib. mass-ft./Ib. force-sec. sq. 
enthalpy, Btu/lb. mass 

heat capacity ratio, C,/C, 

characteristic dimension, ft. 

exponent, defined by Equation (3) 

molecular weight, 1b. mass/Ib. mole 

rotational speed, rpm. 

Eckert number, defined by Equation (10) 

Grashof number, defined by Equation (7) 

Mach number, defined by Equation (12) 

Prandtl number, defined by Equation (9) 
Reynolds number, defined by Equation (8) 
pressure, psia 

volumetric flow, cmf 

gas constant, equal to 1.987 Btu/lb. mole-°R, or 10.73 
psia-cu, ft./Ib. mole-°R 

entropy, Btu/Ib. mass-°R 

temperature, °R 

characteristic velocity, ft./sec. 

velocity of sound, ft./sec. 

polytropic head, ft.-Ib. force/Ib. mass 
compressibility factor 


Greek Symbols 


coefficient of thermal expansion, vol fract/°R 
polytropic efficiency, defined by Equation (4) 
thermal conductivity, Btu/sec.-ft-°R 

viscosity, Ib. mass/ft.-sec. 

kinematic viscosity, sq. ft./sec. 

constant, equal to 3.14159 

density, Ib. mass/cu. ft. 

volumetric flow coefficient, defined by Equation (13) 
polytropic head coefficient, defined by Equation (5) 


Pre ons 


be ee 


Subscripts 
refers to compression stage inlet 
refers to compression stage outlet 
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